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FOREWORD 


This report covers the results of one portion of a two year research 
project carried out by the School of Mechanical Engineering at the Georgia 
Institute of Technology, Atlanta, Georgia for the NASA Manned Spacecraft 
Center, Houston, Texas, under Contract No. NAS 9-10415. This report summarizes 
the results of the radiator simulation analysis. The users manual for the 
execution of the computer program resulting from the analysis in this report 
is contained in a separate report. A third report covers the development of 
a simplified system simulation and the initial phases of a system optimization 
procedure. The project title is "Study of Design Parameters of Space Base and 
Space Shuttle Heat Rejection Systems." The work reported here was monitored 
by Dr. W. E. Simon of the Power Generation Branch of NASA MSC, Houston, Texas, 
and was carried out by Dr. W. Z. Black and Dr. W. Wulff as Co-Investigators 
and Mr. S. M. Morcos, Mr. S. L. Yao, and Mr. R. M. Hinson graduate students 
in the School of Mechanical Engineering under the direction of Dr. S. P. 

Kezios . 

The work carried out by Dr. W. Z. Black is reflected in the analysis 
of Part II, Chapters 2, 7, 8, 11 and 14. Dr. Wulff is responsible for the 
analysis of Part II, Chapters 1 and 3 through 6, Chapters 9, 10, 12 and 13 
as well as the numerical analysis in Part III. Mr. Morcos completed the 
property analysis reported in the Appendices A and B. Mr. Yao completed 
the surface coating property and the shape factor analysis reported in 
Appendices C and D, respectively. 
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SUMMARY 


A transient heat transfer analysis was carried out on a space 
radiator heat rejection system exposed to an arbitrarily prescribed combination 
of aerodynamic heating, solar, albedo and planetary irradiation. A rigorous 
analysis was carried out for the radiation panel and tubes lying in one plane 
and an approximate analysis was used to extend the rigorous analysis to the case 
of a curved panel. For the rigorous analysis the radiator system consists of 
equally spaced parallel coolant flow channels, all in one plane and connected 
by plane fin panels of trapezoidal cross-section, symmetric with respect to 
two normal planes, one passing through the tube axis, the other through the 
center between adjacent tubes. Investigated was one typical tube-fin element 
and the result extended over the entire system, on the basis of the above 
symmetry. 

The rigorous analysis was extended, by approximate methods, to include 
radiator system which do not conform to the above symmetry restrictions. 

As a result, radiator systems can be treated whose coolant flow channels 
do not lie in one plane, provided that the radiative interaction between 
neighboring tube-fin elements is small when compared with the radiant flux 
densitites at the panel. Moreover, radiator systems with non-uniform irradiation, 
non-uniform coolant inlet conditions and U-shaped coolant channels can be 
treated, provided that the spacing between channels is small when compared with 
the channel length. 

The analysis permits the consideration of both gaseous and liquid 
coolant fluids, including liquid metals, under prescribed, time-dependent inlet 
conditions. The flow channels are covered with the same passive thermal control 
coating with optically diffuse but wavelength and temperature dependent 
optical properties. 

The major results of the analysis are the prediction of both transient 
and steady-state, two-dimensional temperature profiles, the local and total 
heat rejection rates, the coolant flow pressure drop in the flow channel, and the 
total system weight and the protection layer thickness. 

A computer program consisting of 62 program units was coded to 
execute the numerical solution of the system of differential equations 


iii 



occurring in the analysis and to predict principal design parameters. The 
modular program structure readily permits later modifications. A separate 
final report entitled "Space Radiator Simulation, Manual for Computer 
Program" has been prepared which describes the computer programs [29] . 

A simplified analysis was carried out to aid the detailed analysis 
and to serve as the basis of systematic optimization. This analysis is 
covered in a companion report, entitled "Simplified Analysis and Optimization 
of Space Base and Space Shuttle Heat Rejection Systems" [30]. Regarding the 
heat rejection rates, its results are, for the test cases carried out so far, 
within 4% in agreement with the results of the detailed analysis. However, 
this rigorous analysis has greater applicability and detail. 

This engineering analysis report is one of three final reports. 

The other two reports are a user’s manual describing the computer code of 
this extensive, rigorous analysis [29] and the final report covering the 
simplified radiator system analyses and system optimization which describe 
both analysis and computer codes [30] . 

This report is written in two principal parts. The analysis and 
the governing equations are contained in the first part, Chapter II, titled 
Analysis. The numerical techniques are covered in Chapter III. Details 
which the reader may need to expand the program are placed in the 
appendices . 
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I . OBJECTIVE 


The purpose of the analysis presented here is to develop a radiator 
system simulation which serves (i) to provide the design parameters necessary 
for the development of the radiator system, (ii) to predict the transient 
radiator performance under prescribed environmental and operation conditions, 

(iii) to predict the system response to conditions which lead to coolant 
fluid temperatures outside their operational temperature range, and (iv) to use 
the system performance data to suggest design options for shuttle radiator 
panel. 

The class of system analyzed here is described in the following 
section. The system parameters produced for design purposes are those 
which describe coolant fluid flow field and the thermal state of the radiator 
structures as well as the geometry and the weight of the system components. 

The arbitrarily prescribed environmental conditions consist of the 
specification of ascent and reentry profiles and of solar, planetary and albedo 
irradiation as a function of time. 

The prescribed operating conditions are the coolant fluid inlet 
properties. 

The analysis is designed to accommodate spectral characteristics 
of surface coatings, specified as functions of temperature. The analysis is, 
however, restricted to optically diffuse coatings. 

Both coolant and structural material properties are accepted as 
prescribed functions of temperature. Coolant properties are specified as 
functions of pressure or density as well as temperature. 

The analysis serves as the basis for a large-scale computer code 
in modular form. Material properties, complex mathematical operations and 
readily identifiable tasks in the computer code are written as subprograms. 

The computer code is designed to simulate space radiator heat rejection 
systems during ascent, reentry and mission phases of the spacecraft and to 
optimize the radiator system configuration via enumeration of parameter sets. 
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The Analysis is covered in the following chapter. The numerical 
techniques employed are discussed in Chapter III, while the preparation for 
the program units describing material properties is deferred to the 
Appendices . 

The computer code is described in a separate manual [29]. 

When the contract began over two years ago, the emphasis on radiator 
design was primarily on the space base heat rejection system. Since that time, 
the emphasis has shifted so that it is now heavily placed on the shuttle 
vehicle. Furthermore, the responsibility of developing the heat rejection 
system has been shifted from the Tower Generation Branch and the design philosophy 
has changed to an integrated system which includes waste heat from sources 
other than the power generation system which was anticipated as the sole 
source of waste heat when the contract began. Due to the shift in design 
philosophy, the report will not recommend detailed design considerations, 
although once heat loads from other sources are known, the analysis presented 
here can be used to aid in the design of an integrated radiator system. 
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II. ANALYSIS 

The analysis is presented in three major parts. In the first part 
are developed the governing equations of transient heat transfer within and 
from the radiator system; these equations are the basis of the numerical 
simulation. The second part of the analysis is devoted to the computation 
of design parameters dictated by operational conditions, while third 
part covers the development of thermodynamic and transport proper- 
ties of structural materials, coolant fluids, and the atmosphere. 
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A. System Description 

Two radiator systems are considered in this analysis. The 
first system considers a flat radiator panel divided by regularly spaced 
coolant channels all having identical inlet coolant fluid properties. 
This system is treated rigorously. The second system considers a 
non- symmetrical radiator panel. The non- symmetrical conditions can be 
caused by a curvature in the radiator panel or by coolant channels 
that are parallel but not equally spaced or formed in the shape of a U. 
This system is treated in an approximate manner. The details of the 
non- symmetrical analysis are given in Section II D. 

For the purposes of system simulation, the radiator system 
consists of four components: 

(i) the fin 

(ii) the coolant fluid 

(iii) the coolant flow channel 

(iv) the meteoroid protection layer 

Two coordinate systems were introduced (see Figure 1), one for the fin 
and the other for the flow channel and the meteoroid protection layer. 
The rectangular Cartesian system (x, z) for the fin has its z-axis 
parallel to the tube axis, starting at the inlet plane, and its x-axis 
passing through the line of profile symmetry, with x = 0 designating the 
root of the fin. Cylindrical coordinates (r,z) are used for both the 
tube and the meteoroid protection layer, with the z-axis along the 
tube . 
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The radiator system is describable in terms of the following 
six dependent state variables : 

(i) for the fin: 

the fin temperature T^(t;x,z) 

(ii) for the coolant fluid: 

the fluid temperature T c (t;z) 
the fluid pressure p(t;z) 

the fluid velocity w(t;z) 

(iii) for the fluid flow channel: 

the tube wall temperature T^(t;r,z) 

(iv) for the meteoroid protection layer 

the protection layer temperature T^(t;r,z) 

These six dependent variables must satisfy four equations of energy 
conservation, one for each component of the system, and further, the 
equations of mass conservation and of momentum balance for the fluid. 

These conservation equations take on the form of partial differential 
equations subject to initial and boundary conditions. Finally, the 
energy conservation equation for the fin involves the net radiant and 
convective power fluxes leaving the fin. 

In the following Sections, II B.2 through II B.5 are discussed, 
in that order, the six principal governing equations associated with the 
four system components listed above, the subsidiary equations governing 
the radiative heat exchange and, lastly, the convective heat transfer 
between coolant fluid and tube wall and between the fin and the atmosphere 
during ascent and reentry. 
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B . HEAT TRANSFER 

1 . Introduction 

The objective of the radiator simulation is to predict both the 
transient and the steady state heat transfer characteristics under pre- 
determined operating conditions, stationary in the latter case and dynamic 
in the former. Both cases were treated as initial value problems, and the 
principal governing equations are partial differential equations which are 
linear in the time-derivatives of first order at most. 

Under stationary boundary conditions steady state will be reached, 
regardless of initial conditions*, as all partial derivatives with respect 
to time vanish on account of dissipative effects within the system. For 
the computer simulation of steady state conditions this means that the 
process of advancing in time can be discontinued as soon as all variables 
y., i = 1, 2, . . . N, have reached their expected asymptotic values (y^)^ 
with sufficient accuracy, that is when for some chosen £i 



has been reached. The expected asymptotic values (y can be estimated 
on the basis of the recognition that for large enough values of the time, 

6 > "m 

y i ■* ± ae bt 0-2) 

with, a > 0, b > 0, to be determined; t > The evaluation 6 ^ during 

the numerical integration is covered in Chapter III. 


*Subject to certain continuity requirements which are discussed in Chapter III. 
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2. The Fin 

The objective of this section is to develop the energy equation 
for the fin. The derivation is based on the assumptions that the thermal 
conductivity and specific heat of the fin material are functions of 
temperature while the fin density is constant. The energy balance for 
the fin accounts for both radiative and convective fluxes from the fin 
surfaces. The development of a method to predict the net radiant flux 
from the fin surfaces can be found in Section II-6. The procedure used 
for the evaluation of the convective flux from the fin surface is outlined 
in Section II-7. 

The energy balance on a differential volume of the fin can be 
expressed as 


3_ 

3x 


3T 


k„A 


f x 3x 


dx + 


3z 


k A 
f z 


ST^ 

3z 


dz + q’ 


+ 


net , rad 


ti 

q A 

aero s 


p f Vc P f 



( 2 . 1 ) 


where T is the fin temperature and the coordinate system is shown in 

Fig. 1 . The areas A and A represent the cross-sectional areas of t 

x z 

fin perpendicular to the heat flow in the x and z directions, and A 


represents the total non- adiabatic surface area of the fin element. 

The symbol V is the volume of the fin element. The properties of the fin 
material are represented by the symbols k^,p^ and c ^ which stand for 
the thermal conductivity, density and specific heat, respectively. The 
terms q" , and q M appearing in Eq. 2.1 denote the radiative and 

I"1 O t j IT did d6 IT o 

convective heat gain from the surroundings to the fin surfaces. 


The first two terms in Eq. 2.1 constitute the net conduction 
of energy into the fin element. The third and fourth terms on the left 
hand side of the equation stand for the net radiation and convection gain, 
respectively, from the surfaces of the fin element. The tem of the right 
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hand side of the equation represents the storage of internal energy within 
the fin element. 

The appropriate areas and the volume to be substituted into 
Eq. 2.1 are 


dA x = 2[s^ - ex] dz 


( 2 . 2 ) 


dA^ = 2[s^ - ex] dx 


(2.3) 


dA g = dxdz / [c + l] 


1/2 


(2.4) 


dV = 2[s^ - exj dxdz 


(2.5) 


where s r is the fin half thickness at its root and c is the negative slope 
f 

of the fin side surfaces. 

Substituting Eqs. 2.2 through 2.5 into Eq. 2.1 and simplifying 
the resulting equation yields 


8T, 


f >f 3 1 ~ k f 


Pr C 


z 


r- 



3 T f + 3^ 

9x 2 3 Z 2 

dk_ 

+ — - 
dT 

(m 

y 3x 

Mi 

ij 



- 




-(&) 


_3T 

8x 


vx 


+ 2(t-cx) ^ q net,rad + q aero^ 


( 2 . 6 ) 
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The Normalization of Eq . 2.6 is achieved by defining nine dimension- 
less quantities: 

Let 


£ = x/H 


(2.7) 


C ~ z/L 


( 2 . 8 ) 


represent the nondimens ional fin coordinate perpendicular and parallel to 
the tube, respectively. Then 


t = tw /L 
o 


(2.9) 


is the dimensionless time. The symbol w 0 stands for the velocity of the 
coolant fluid entering the tube. 

The nondimens ional fin temperature is defined as 


T f (t,C,C) 

e f ( T ,£,0 - 

o 


( 2 . 10 ) 


where T q is the temperature of the coolant fluid entering the tube. 

The dimensionless conduction parameter N^ c and the Fourier 


number are defined as 
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N = HaT J /k, 
Inc of 


( 2 . 11 ) 


N 


a f L _ k f L 


Fo „2 2 

w H p r c r w H 

o f pf o 


( 2 . 12 ) 


and the nondimens ional geometrical quantities are defined as 


S f = s f^ H 


(2.13) 


B(4) = 2 (s f - cO 


D = 


(c 2 + 1) 


1/2 


(2.14) 


(2.15) 


where c is the negative slope of the fin side surfaces. For a non-tapered 
fin c is zero. 

Both the net radiative and convective flux terms appearing in 

4 

Eq. 2.6 are normalized by dividing each term by aT , or 


*Vet » ^ad = 


Siet ,rad 
4 

qT 


l aero 


aero 


aT 


(2.16) 


(2.17) 
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The energy equation, Eq . 2.6, may now be written in terms of 
the nondimensional quantities given in Eqs. 2.7 through 2.17. The resulting 
nondimens ional equation is 


8 f “ N Fo < + 




(fj 


T dk, 

(e f ) +^7— 

f c? k f dt 


(0^) + 


(?/ 


(e f ) 


C J 


(2.18) 


B '”f' £ ■ B "Nelnet, rad + W,conv ] j 


- 2 I (8.) +7 N„Jq 


The dot superscript appearing in Eq. 2.18 denotes differentiation 
with respect to nondimensional time and the subscripts £ and £ represent 
partial differentiation with respect to the dimensionless coordinates 
indicated. 

The normalized energy equation, Eq . 2.18, defines the rate of 
change of the dimensionless fin temperature 0^. 

The Boundary Conditions for Eq. 2.18 are taken as follows: 

(i) The fin root is at the temperature of the outside of the 
tube . 

(ii) The fin tip is insulated. 

(iii) The portion of the fin in contact with the inlet manifold 
is at the outside temperature of the manifold. 

(iv) The portion of the fin in contact with the exit manifold 
is at the temperature of the outlet manifold. 

Written in mathematical terms these four boundary conditions are 

G f (T,0,c) = 9 w (t,£ o ,£) (2.19) 





(t, 1,£) = o 


( 2 . 20 ) 
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e f (T,5,o) = e w (T,5 o ,o) 


( 2 . 21 ) 


e f (T,5.i) - e (t,c ,i) 

■L WO 


( 2 . 22 ) 


Implied in Eqs. 2.21 and 22 is, firstly, that the fluid temperature 
in the manifolds be equal to the fluid temperature at the tube inlet (inlet 
manifold) and at the tube exit (exit manifold) and that it remain unaltered 
along the manifold and vary only in time; secondly, that the temperature 
drop through the manifold wall be equal to that through the tube wall. 

The latter assumption is well justified because the temperature drop is 
exceedingly small. 

It may be noted in connection with the boundary conditions that 
the radiative interactions between manifold and fin as well as between 
manifold and tube are not taken into consideration at this time. 

The Initial Condition for Eq. 2.18 may be any arbitrary relation 
representing a continuous temperature distribution over the fin, including 
the boundaries. The selection of the initial fin temperature distribution 
is left to the user. 
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3 . The Coolant Fluid 

The objective is to develop a unified treatment of all possible 
coolant fluids, that is gases, dielectric fluids and liquid metals. Three 
principal governing equations are sought which, together with the necessary 
thermodynamic and transport properties specified for each fluid of interest, 
define the fluid temperature T , the fluid pressure p (and thus the thermo- 
static state of the fluid) and the fluid velocity w, all as functions of 
time t and axial distance z. 

The Continuity Equation for one dimensional flow through channels of con- 
stant cross-sections reads 

+ TT (P w ) = 0 (3.1) 

d t d Z 


where p represents the fluid density. Replacing the density through the 
thermal equation of state 


P = P (P> T c > 


(3.2) 


renders the continuity equation in terms of derivatives of the primary 

variables T , p and w: 
c r 


3 



3w 

3z 


+ w 


’(* 


3Tc 

3z 


-k 



(3.3) 


where k and 6 stand for the isothermal compressibility and the isobaric 
expansion coefficient, respectively 


K 



(3.4) 
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(3.5) 


The Momentum Equation 






+ pa z 


(3.6) 


constitutes the balance between inertia forces on the left-hand side, 
pressure forces, wall friction and external field forces (gravity) in axial 
direction, on the right-hand side of Eq. 3.6. The Fanning friction factor 
f is a function of the Reynolds -number = d w/v, where d and v stand 

for the tube diameter and the kinematic viscosity, respectively and sub- 
script o designates the fluid inlet conditions. The following relations 
are used to compute the Fanning friction factor: 

for < 2300 (laminar flow) 

Re - - 

4f -4*- (3-7) 

N Re 


for 2300 < N < 10 6 
“ Re — 


(Ref. 1) 


4f 


0.0054 + 0.396 N 


- 0 . 

Re 


3 


(3.8) 


for N Re > 


10 (Ref. 2) 


4 f B 0.0032 + 0.221 


N 


-0.237 

Re 


(3.9) 


Equations 3.8 and 3.9 could be replaced by a single equation (Ref. 3) 
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4f = [0.86859 Jin 


Re 


1.964 £n N 


Re 


3.8215] 


(3.10) 


which, however, requires more computational effort. 

The Energy Equation . Let u, h, h^ and = T^(t; r_^,k) represent, 
respectively the internal energy and the enthalpy of the fluid, the convec- 
tive film coefficient and the tube wall temperature at the fluid-wall 
interface. For one-dimensional flow through channels of constant cross- 
section and with heating or cooling from the wall, conservation of energy 
requires that, with respect to a stationary reference frame 


^ tP(u + 4] [p« < h+ 4 


- a z) ] = 7 h (T - T ) 
z d c w c 


(3. IX) 


+ (k 

d Z C 


8T 
( 

dz 


The first term constitutes the storage of thermal and kinetic energies; 
the storage of potential energy, being negligibly small for expected 
accelerations, is ignored. The second term on the left-hand side stands 
for the convection of thermal, kinetic and potential energies as well as 
the power associated with the pressure. The right-hand side contains, 
firstly, the convective heat transfer from the channel wall to the fluid 
and, secondly, the axial heat conduction term which will later be shown 
to be negligibly small for all fluids. The factor 4/d in front of the 
convective heat transfer term results from the ratio of the channel 
circumference to the cross-sectional area, evaluated for a circular tube. 
The symbol k c represents the thermal conductivity of the coolant. 

Given a caloric equation of state 


u = u(p,T ) 
c 


(3.12) 
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or an equivalent expression for the specific heat, c^, one can write, 
with the aid of Eqs. 3.2, 4 and 5 


dh = c dT 
P c 


+ J (1 - gT c ) dp 


(3.13) 


du 


(c - dT + — (icp - $T) dp 
P p c p ^ J v 


(3.14) 


and then recast Eq. 3.11 in terms of the derivatives of the principal 

variables T , p and w: 
c 


P c 


3T 

c 

9t 


- 6T 


c 9t 


+ pw 


3w 

3t 


f h 

d c 


(T - 
w 


V 


+ 


3z 


(k 


3T 

c 

3z 


) - pw 




+ 



f3T) 


iE 

3z 


+ w 


3w 

3z 


” a 

z 


] 


(3.15) 

Thus, three governing equations, Eqs. 3.3, 3.6 and 3.15, have been 
established which define the three principal variables T , p and w, 
provided the initial and boundary conditions are properly specified and 
the convective film coefficient can be predicted. Since the normalization 
of these equations, followed by an order-of-magnitude comparison will 
indicate that the conductive term in Eq. 3.15 is negligible so as to 
simplify the boundary conditions, the discussion of the initial and 
boundary conditions is deferred until after the normalization. 

The convective film coefficient is computed from the following 
relationships between the Nusselt number = h^d/k^, the Reynolds number 
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N and the Prandtl number 1SL = u c /k where y represents the dynamic 
Re Pr p c J 

viscosity: 


For N_ <0.1 (liquid metals) 

Pr 

N Nu = 6.5 + 0.025 (N Re N pr ) 0 ' 8 


(3.16) 


which produces Nusselt numbers between those appropriate for uniform 
heat flux (Martinelli) and for uniform wall temperature (Seban and 
Shimazaki) (Ref. 4) . 


For N_. >0.1 and 

Pr 

for N_, < 2300 (laminar flow. Ref. 5) 

Re 

0.0668 (N Re N pr 6) 

N Nu = 3-65 + 1 + 0.045 (N Re N pr «) 


(3.17) 


for KL 2 l 2300 (turbulent flow, Ref. 5) 
Re 


N = 0.116 (N 


0.667 _ 125) N 0.333 


(1 +«) 


(3.18) 


where 6 = d/L stands for the tube diameter-over-length ratio. 

The Normalization of Eqs. 3.3, 6 and 15 is carried out for the purpose 
of scaling and performing an order-of-magnitude comparison. The 
computational effort is also reduced in the process. 

Let 


C = z/L 

= tw /L 
o 


(3.19) 


T 


(3.20) 
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represent the nondimens ional axial distance and the nondimens ional time, 
respectively, and let the dimensionless state variables be defined as 


0 C <T,S) 


T c (t ’ z) 

T 


(3.21) 


7T (T,C 




(3.22) 


on (t,<; 


W (t.z) 
w 

o 


(3.23) 


and to represent the nondimens ional temperature, pressure and velocity 
of the fluid, that is, the principal dependent fluid flow variables. 

The subscript "o" designates the constant reference state of the fluid 
at the tube entrance. Introduce next the nondimens ional density 

v( T , £) = --P ^ * z ^~ (3.24) 

p o ’ 

the nondimensional radial distance from the channel axis 

2r 

n - (3.25) 

and the nondimensional tube wall temperature 

T (t;r,z) 

e (i;n,c) = (3 * 26) 

w I 

o 

Notice, that all dependent variables lie between zero and unity, except 
v and w whose product (vw) remains essentially equal to unity with neither 
v nor oa departing far from unity. The reference temperature T q is, under 
normal operating conditions, the highest temperature in the system. 
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Furthermore, consider the following (|>-values to vary along the 
channel axis: 


= T 


= p K 


i C 

<P = £ 

cp c (p T ) 
p o o 


c 

„£. 


c 

p,o 


k = 

k (p T ) k 
o o o 


and finally, the following constant F-parameters : 


F 

P 



F 


f 



4f 

6 


F 

z 



P 


c T 
° p,o o 


2 

w o 

F = F /F = 

M z'p c T 
Po o 


(3.27) 


(3.28) 


(3.29) 


(3.30) 


(3.31) 


(3.32) 


(3.33) 


(3.34) 
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Let the dot above a variable designate partial differential with respect 
to the nondimens ional time t and the subscript £ partial differentiation 
with respect to the dimensionless axial coordinate £ . Then the principal 
conservation equations, Eqs. 3.3, 3.6 and 3.15 read, respectively and in 
nondimens ional form: 


„ [0 + a)(e ) 

3 c c £ 


w + d> (tt + am ) 
C k C 


(3.35) 


0) -{- 0) 0J^ 



(3.36) 


0 


4N 


Nu 


®c - “V% e N pr • v* cp ~c' ' 4 N Nu 


1 (0 - 0 ) + — - 


cp 


cp 


*u K <V, \ -“ ( V C + 


1-0 <h . , 

r — ^ ~ it + f m -f-"" 

Z *cp V ? M ♦cp 5 


(3.37) 


These are the three equations which define the three time rates of change ,l 

• • • 

0 , 7 T and m. It can be seen' from Eq. 3.37 that the axial conduction 
c 2 

is always small, of an order less than 6 (since > 3) s when compared with 

the convective term, unless the fluid should reach the wall temperature 

within the very first small fraction of the tube length. Axial conduction 

r\ /" 

is hence ignored as 6 * 10 , and the order of differentiation of Eq. 3.37 

is reduced to one. 
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. d> 0 „ 4N 

q _ p — & £. ' j -p 1 • _ Nu 

c z (p v M $ wt0 6N N 

C P r cp Rc Pr 


0-0 

w c 


V 0 


cp 


- U',{(9 ) + F x 

c c z 


1-00 F 

£& I Ji WM 

, TT * r 

vd) c Y 

r cp ^ cp 


(3.38) 


The Boundary Conditions to be imposed on Eqs. 3.35, 3.36 and 3.38 are chosen, 
at the channel entrance, to be 

(i) mass flow rate m, prescribed as a function of time 

(ii) constant inlet pressure P q . 

(iii) continuous transition, of the inlet fluid temperature, from an 

initial temperature 0 = 0 ^ to the constant operational temperature 
0(t,O) = 1, or an arbitrarily prescribed inlet fluid temperature that 
is a continuous function of time 

These boundary conditions accommodate the calculation of the steady state 
conditions as well as of the most likely start-up operation toward stationary 
operating conditions. Notice that there are no step changes implied in 
any of the dependent variables which is essential for the numerical inte- 
gration. Writing these boundary conditions more specif ically, one gets 
at 4 = 0 

9 (t, 0 ) = 1 - (1 - 6 .) e~ ?T ( 3 . 39 ) 


7T (l, 0) = 1 


(3.40) 


W (T, 0) 


1 

v(tt, 0) 


(3.41) 


Here, the time constant was chosen arbitrarily so as to have the fluid 
reach, within 0 . 1 %, its steady inlet temperature at the inlet during the 
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time interval that is takes a fluid particle to pass through the channel. 
Equation 3.41 is given through the thermal equation of state, Eq. 3.2. 


The Initial Conditions appropriate to the system of Eqs. 3.35, 3.36 and 

3.38 are derived from the requirement that the flow should initially be at 

steady state with the fluid inlet temperature equal to the uniform channel 

wall temperature. Setting the time derivatives equal to zero in Eqs. 3.35, 

3.36 and 3.38 results in a system of 'three*' ordinary differential equations 

which are linear in dir/d^, d0 /d£ and dm/d^: 

c 


F 

P 


dji 

dC 


dco 

d? 



(3.42) 


^8 dn 
F — ^ j — 
z v dc 


+ 


d 0 
c 

d? 



w 


doo 

d? 


4N. 


Nu 


<5N_ N 
Rc Pr 


(0 “ 9 > 
w c 


(3.43) 


, dO 

, dTT a C 

d 5 \ dc 


1 du> 
w d£ 


= 0 


(3.44) 


These equations can be solved subject to 


_ 1 

w v(n, 6 ) 
c 

provided the function 0 =0 (0; 1,£) prescribing the initial channel 

r w w 

wall temperature is specified. For the present analysis 0^ was set equal 
to 6.. Equations 3.42 through 3.45 define the initial flow field. 

Quasi- Steady Flow . It may be recognized that the momentum transport 
takes place at a much smaller time scale than the transport of thermal 
energy in that the pressure and the velocity fields adjust virtually 


the boundary conditions at 


(3.45) 
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instantaneously to a change in flow inlet conditions while the response 
of the temperature field to a change in channel wall temperature takes 
considerably longer, the reason being that the pressure perturbatxons 
propagate along the channel with the speed of sound. Unless one is 
specifically interested in the motion of sound waves one may consider 
the dynamics of the flow field, that is the pressure and velocity distri 
butions, as part of the boundary conditions and imposed instantly and 
adiabatically by the flow inlet conditions. 

The fluid temperature remains stationary during the dynamic 
adjustment and the time rate of change of both pressure and temperature 
remain small since ordinarily the pressure gradient remains balanced by 
the wall shear (and by the convective acceleration in the case of a 
gaseous coolant medium). Consequently, Eqs. 3.42 and 3.44 may serve to 
establish the pressure and velocity fields at all times, subject to 
boundary conditions given by Eqs. 3.40 and 3.41 while the temperature 
field remains defined by Eqs. 3.35, 36, 38 and the initial conditions 
discussed above. Particularly, solving Eqs. 3.35, 36 and 38 for 0 c 
gives the differential equation which governs the temperature field. 


0 = 
c 


A Vi 

<b <t> v 
r cp K 


1 - 


4N 


Nu 


6 N p N p 

Re Pr 


(0 


w 


- 0 ) 
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+ F 


Xp 


x [ j> e 


-”«v c + Vf 


cp 


(3.46) 
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This completes the discussion of the development of the governing 
differential equations for the coolant fluid. All equations are solved 
numerically as discussed in Chapter III. The thermodynamic properties c^, 

3 , and « are derived from the thermal equation of state, Eq. 3.2 and from 
the zero-pressure specific heat or other available properties. All thermo- 
dynamic functions as well as the transport properties k and P are considered, 
in general, as functions of two state variables, that is, of p and T or of 
p and T, as discussed in Section E of Chapter II. 
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4. The Flow Channel 


The flow channel is treated as a circular tube with inner radius 
r^ = d/2 and outer radius r Q = r\ + s^. The tube wall temperature 
T w (t;r,z) is defined through the familiar equation of energy conservation, 
written for the case of circular symmetry: 


3T 

w 

at 


1 a 


a r r 




3 2 T 


“w [ r If (r i^ } + 


w 


] + 


az 




(4.1) 


where a represents the thermal diffusivity of the wall material and 
w 

(pc) its volumetric heat capacity. The boundary conditions are 
w 

at r = r. 
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h 
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(T - 
w 


T )= k 
c' w 


3T 

w 

ar 


(4.2) 


at r = r 
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3T 


w 


w 3r 


- rpk. 


3T^ 

lx 


(2tt - ij j) k 
T m 


aT 
m 

ar 


0 


(4.3) 


Here 

h is the convective film coefficient, 
c 

T c the fluid temperature, 

k the thermal conductivity, 

ijj the portion of outer tube circumference in contact with fins, 

expressed in radians, 

and the subscripts w, f and m designate, respectively, tube wall, fin and 
meteoroid protection layer. Equations 4.2 and 3 constitute the continuity 
of the heat flux at the fluid-wall interface and at the wall-fin and wall- 
protection layer interfaces. Circumferential temperature variations are 
ignored. 
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After introducing 


n 




dk a L 

w _ w 

dT » ^Fo ,w 2 

w r . 
o 1 



(4.4) 


with H representing the fin height, that is, the distance between the fin 
root and fin tip, one may recast Eqs. 4.1, 2 and 3 to read 


®w ^Fo,w 


1 

n 



(V 
^2 -- 


8 2 G 


w 





(4.5) 


(4.6) 



28 


at n = n 


o 


30 30 , 

w r 

— - w 


30 

(2tt - ip) <J) = 0 

T m 3n 


(4.7) 


The superscript dot represents partial differentiation with 
respect to nondimens ional time t . As before in the treatment of the coolant 
fluid one recognizes that, with 6^-10 , axial conduction remains 

insignificant. Thus, Eq . 4.5 takes on this final form 


0 

w 


Fo ,w 


1 

n Sn 


30 


(n 
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dX] 


) + 



(- 


n 

3n 



(4.8) 


Equations 4.6, 7 and 8 define the tube wall temperature, provided Eqs. 4.6 
and 7 hold at some initial time and the initial temperature 6 w (0;n,C) is 
prescribed as a sufficiently smooth function of n and £ . 


Low Biot Number . When one compares possible heat fluxes at the 
fluid-wall interface with the possible radiant fluxes from the outer 
surface of the meteoroid protection layer covering the tube, one con- 
cludes that the maximum fluxes occur at the inner tube wall and that the 

Biot number . in Eq. 4.6 is the largest ratio to be expected of external 
Bi 

to internal thermal conductances. Thus, if N.^ is small, say less than 
0.05 (Ref. 6), then the temperature variation inside the tube wall is too 
small for experimental detection and a computation of the detailed 
temperature distribution on the basis of Eq. 4.8 cannot be justified 
as the associated computational effort is considerable. 

In cases where the equivalent Biot number, representing the 
total thermal resistance within the channel wall and the protection layer, 
is small, the tube, the protection layer and a representative portion of 
the fin root are combined into a single control volume as depicted in 
Fig. 2. The equivalent Biot number N and the chosen limit are 


N Bi N Nu d 


s s 

(r 2 + 7 s ) < 0.05 
k k 

w m 


(4.9) 
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where s and s stand for, respectively, the tube wall and the protection 
w m 

layer thicknesses. The combined volumetric heat capacity per unit of 
axial distance for the control volume consists of three parts, the first 
one for the tube wall, the second for the protection layer and the third 
for the first half fin element: 


( s 

nd (pcs)^ j (1 + ~~) + (1 + 2 
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HU 


(pcs) 


d + d ^ (pcs) 


m Ax 

-rrd 

w 


w 



s (pc) 

r )] -&T 

r T w 


= X_ Trd (pcs) 
f ' w 


(4.10) 


where the subscripts w, m, f, r and t designate, respectively, parameters 
of the tube wall, the meteoroid protection layer, the fin, the fin root 
and the fin tip and where 

d is the tube diameter, 

p the density, 

s the thickness, 

Ax the node spacing on the fin, and 

A£ = Ax/H, the nondimens ional node spacing 

Heat enters the control volume, per units of time and axial 
distance, by convection from the fluid 

Trd h (T - T ) , 
c w c 

and by convection and/or radiation from outside 

(2tt - (f + s w + sj + Ax (^ u x + q M 2 ) f 
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where the subscripts 1 and 2 designate upper and lower fin sides, 
respectively, and where q ,f represents the sum of the convective and of 
the net radiant fluxes entering the outer surfaces. In the case of q^, 
the average over upper and lower portions of the outer circumference is 
to be taken. 

Heat leaves the control volume through the fin by conduction, 
again per units of time and axial distance 


2 (s 


8T* 

3x ; f 


evaluated at x - Ax/ 2. 

Combining the last four expressions into the energy balance 
leads to this expression 
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(s 
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(4.11) 


where is defined by Eq. 4.10 and 


1 )J. .1 S W + s m. 

" (2 - P ( I + ~T^ ) 


( 4 . 12 ) 


Equation 4.11 replaces Eqs. 4.6, 7 and 8 as well as Eqs. 5.3, and 5 
governing the temperature distribution in the meteroid protection 
layer and, lastly, Eq. 2.19 which constitutes one of the boundary 
conditions for the differential equation governing the fin temperature 
distribution. The only condition under which all the above equations 
may be replaced by the single equation, Eq. 4.11, is given by Eq. 4.9. 
Finally, it may be noted that Eq. 4.11 could also be normalized but since 
no new dimensionless groups result from such normalization it is omitted 


here. 
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5 . The Meteoroid Protection Layer 

The differential equation governing the temperature distribution 
within the meteoroid protection layer which covers the flow channel, is 
identical to that for the channel wall, Eq . 4.8, except for the two 
dimensionless parameters, the Fourier coefficient <}) and the conductivity 
temperature coefficient cj)^ which now must be evaluated for the protection 
layer material: 

a 

d> = — <j> 

Fo,m a w Fo,w ( 5 . 1 ) 



dk 
m 

dT 


(5.2) 


After ignoring the axial conduction for the reasons stated in Section 4 
one obtains as the nondimens ional energy conservation equation 


0 


m 





(5.3) 


Two boundary conditions are required, one of which is given by Eq. 4.7 
while the other one is dictated by the heat flux continuity at the outer 
boundary : 

at r = r 
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BT 


m 


dr 


m 


(5.4) 


where q" is defined as the net flux entering both by radiation and/or 
aerodynamic heating. Equation 5.4 reads in nondimensional form at n = 


— 2- =4) (q + q ) 

3n Nc M net,rad aero 


(5.5) 
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with 4 > Nc representing a local conductance parameter 
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odT 


Nc 2k 
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(5.6) 


(5.7) 


(5.8) 


Here, a represents the Stefan- Boltzmann constant, and q" , and q" 

110 L y 1. 3 . 0 . W0i o 

the net radiant incident heat flux and the convective heat flux, respectively. 

If the initial temperature distribution is given, and if Eqs. 4.7 
and 5.5 hold initially, then the protection layer temperature 6^ is completely 
defined as a function of time and location by Eqs. 5.3, 5 and 4.7, provided 
the incident radiation and the aerodynamic heating are prescribed as 
functions of time. The prediction of these incident heat fluxes is dis- 
cussed in Section 6 and 7 of this Chapter. 

In the case of low Biot numbers at the coolant fluid-tube wall 
interface the meteoroid protection layer is lumped together with the tube 
wall as discussed in Section 4 of Chapter II. 
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6 . Radiation 

The objective of this section is to develop a procedure 
to predict the net radiant heat flux incident on fin and tube surfaces 
which are exposed to any combination so solar, albedo and planetary 
irradiation. Included into the assessment of radiative transfer is the 
radiative energy exchange between fin panel and flow channel or its 
protective coating as well as the effect of structural panels in the 
vicinity of the fin system; however, not included is any possible gas 
radiation as could conceivably be encountered during reentry. 

In seeking the proper mathematical model it is recognized, 
firstly, that the prevailing thermal radiant energy lies in either the 
visible (solar irradiation) or the infrared portion of the spectrum, and, 
secondly, that the fin system is coated, for the purpose of optical 
optimization, with a dielectric paint. Consequently, spectrally dependent 
optical properties must be dealt with, but, for the latter reason, the 
transfer matrix of radiative exchange can be expected to remain temperature 
insensitive over some range of operational conditions, a fact which is 
very much appreciated from the computational view point. 


For the analysis, the fin surface A^, the outer surface on the 

flow channel A , and the structural surface (s) A , are all considered 
n n 

as parts of an enclosure C which is completed by a set of arbitrarily 

concave, non-reflecting imaginary surfaces A^ which connect A^, A^ and 

A^ and along which is specified the emerging net radiant heat flux 

representing solar, albedo and planetary irradiation. The sum A^ + A^ + 

A + A is the inner surface A of the enclosure, 
n e c 

Three steps are necessary for the prediction of the incident 

net radiant heat flux q M .. Firstly, the elemental exchange areas* 

^net,rad 


fi!i = 


3A. 3A. 
i J 


CQS<Pj- COS 9 j 


irr 


( 6 . 1 ) 


*For terminology and notations consult Ref. 7, Chapter 2. 
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need to be computed on the basis of the geometric relation between fin 

panels and flow channels. Here, the symbol r designates the distance 

between two area elements dA. and dA. which are visible from each other, 

i J 

and <J> . and <j> . represent the respective angles between r and the surface 
J 

normals on dA. and dA,. The next step is to compute, from its definition, 
the radiosity or leaving radiant flux density, : 


W. 

J 





dA + 


( 6 . 2 ) 



c 


where W , E and e . , stand for the monochromatic radiosity, the 
j,A j,A 

monochromatic black body emissive power 



2 2 

2 he n 1 

5 (hc)/(kXT ) 

■ e 3 - 1 


(6.3) 


and the monochromatic hemispherical emittance, respectively; the subscripts 
i and j designate two discrete points on A^, and A represents the wave- 
length. The Eq. 6.3 constitutes Planck's law of monochromatic emissive 
power intensity; h, c and k stand for, respectively, Planck's constant, 
the speed of light in vacuo and the Boltzmann constant. The third and 
last step is to calculate, on the basis of local energy balance, the net 
incident heat flux 


^net ,rad^ . 
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iot 1 dA : 
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(6.4) 
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It should be obvious that Eqs. 6.1, 2 and 4 contain all the 
necessary fundamental principles but their evaluation will introduce a 
number of simplifications and modifications, each selected for the particular 
system of interest. Specifically, Eq . 6.4 will have to supplement Eq. 6.2 
for portions of where the heat flux is specified. More importantly, 
however, there is a choice to be made in view of the computational process 
regarding particularly Eq . 6.2. One may either solve the monochromatic 
version of Eq. 6.2 n times for the n significant spectral intervals 
encountered and thus face the ultimate task of solving n x N simultaneous 
linear algebraic equations when N discrete points on need to be con- 
sidered (possibly at several time steps during the calculation process) 
and then integrate the resulting total interchange areas over the spectrum 
(see Ch. 5.6 of Ref. 7). Or, one may force the non-gray surface analysis 
into a gray surface analysis by placing the burden of complexity on the 
evaluation of appropriate optical properties. Both techniques afford 
any desirable accuracy of allowing for the spectral differences in surface 
properties, limited only by available calculation time; but the latter 
technique was chosen because, as a result of this choice, the complexity 
remains at peripheral parts of the computer code which are more accessible 
for later modifications toward greater sophistication, also the complexity 
may turn out, in almost all cases, to reduce partly to simple hand 
calculations . 

After introducing 
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Equation 6.2 simplifies to 
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which reduces to the gray-surface radiosity equation whenever Eq. 6.6 

reduces to ot_ = e^. The difficulty now lies in evaluating Eq. 6.6 even 

though the radiosity W. , is as yet unknown. 

J > A 

By successively substituting the right-hand side of Eq. 6.2, 

in its monochromatic form, for W. , on the right-hand side of that 

i > x ^ 

equation, one obtains first (Ref. 8) 


W. = e E + (1 - e 
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( 6 . 8 ) 
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and subsequently a . . as the quotient of two infinite series. Since the 
enclosure radiation is dominated by the fin-sun and fin-sky interaction 
and since Eq. 6.8 contributes significantly only to the fin-flow channel 
interaction, the infinite series in Eq. 6.8 may be terminated after two 
terms (two reflections; the resulting error is less than the uncertainty 


in e ), and Eq . 6.6 becomes: 


X. . E. + 
ij j A 
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(6.9) 
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where 


X. . = 

ij 




( 6 . 10 ) 


( 6 . 11 ) 


In cases where the net radiant flux is specified over portions of A, the 
emissive power E is to be replaced by the net radiant flux q" in Eqs. 6.9, 
10 and 12, which results in one additional term each in the numerator and 
the denominator of Eq. 6.9. 

In summary, the incident net radiant heat flux for the diffuse, 
non-gray enclosure is calculated on the basis of an approximate gray- 
surface analysis in accordance with Eqs. 6.4, 7 and 9 through 11. The 
spectral differences of the surfaces are accounted for in Eqs. 9 through 
11. The remainder of this section is devoted to the solution of the 
radiosity equation, Eq. 6.7. 

Recalling that A is the sum A + A + A_ + A of the outer 
c m n f e 

channel surface A , the possibly present, nearby structural surfaces 
m 

A , the fin surface A,., and the remainder of the enclosure A , one 
n f e 

recognizes that the integrals over A^ in Eqs. 6.4, 7 and 9 need to be 
evaluated twice for each of the four parts, namely once with j = 1 
representing the fin area and then with j = 2 for the exposed channel 
area. Since the incident solar, albedo and planetary radiant flux 
intensities are uniform over the fin area and averaged over the circum- 
ference of the channel area 


[ ^ s i s l 
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( 6 . 12 ) 
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3 s i s 2 
3A i 3A 2 


e 


(1 - a. 2 ) 


W.dA. 

l l 


(1 


a e2 ) q 


fl 
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(6.13) 


where q" designates incident solar, albedo and planetary heat fluxes, 
appropriately averaged over a chosen area element. Should any structural 
surfaces obstruct the incident radiant fluxes (A^ 4 0) then the right- 
hand sides of Eqs. 6.12 and 13 would have to be modified and reduced in 
the shaded portions of and A^; and, if there are m such surfaces. 
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(1 - « kj ) w k 


n 


(6.14) 


j - 1, 2, ... m + 2. 


Obviously, the radiosity and the temperature are assumed to be uniform 
over each structural component. No such structural components were 
included in the analysis reported here, and Eq. 6.14 is taken to be zero. 


This leaves only the integrals over A and A,. to be discussed. 

m r 

Moreover, since Eq. 6.1 is symmetric with respect to its subscripts i 
and j, the elemental exchange area is to be evaluated only once. 

Considering first the fin, that is j = 1 and i - 2, and the 
fact that over the channel surface the radiosity and the temperature 
are considered to be functions of axial distance only 

f 8 2 S l S i 

I 3 A 3A. ^ " a li^ W ± dA i = 

A 

2 



40 


Z=L <j> 

jll - a li ( z ) J W.(z) j ^ -Li j Rd4> dz = ( 6 . 15 ) 

2=0 <j)=o 

z=L 

/u- a 1:L ( z )] W ± (z) SS(z;X fJ z f ) dz 
z=0 


where L designates the tube length; R and are the polar coordinates of A^,, 
with origin on the tube axis, with <p = 0 and <p = representing, respectively, 
the root of the fin and the contact line between the tube and its tangent 
plane through the center of AA^ on the fin. The first step in Eq. 6.15 was 
obtained by integrating over AA^ and subsequently applying the mean-value 
theorem of integral calculus, while the second step simply defines the ex- 
change function for every point (x^, z^) on the fin which was integrated in 
closed form for the right-circular flow channel. The result is shown in 
Appendix D. 


The exchange function of the tube with respect to the fin is ob- 
tained by dividing SS in Eq. 6.15 by (R^*) . Thus 
,2 


3A„3A. 

2 i 


/ 

A L, H. 

1 *// 


(! - o ) w. dA. = 


[1 - aU,x., 2 .) ] W(x f ,z f ) SS (z ,x., 2 .) dx. dz f 


( 6 . 16 ) 


z=0 x=0 


For the numerical evaluation of the integrals a suitable quadrature such as 
the trapezoidal rule is chosen so as to render Eq. 6.7 in this form 
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P. = Y M. . W. (6.17) 

J j Ji 1 

i,j = 1, 2, . . . , N 


which is a system of N linear algebraic equations for the N = (n x + 1) (n g + 2) 

unknown values of the radiosity W. . Here, n and n are the numbers of sub- 

X X z 

divisions chosen in the x- and the z-directions , respectively. The vector P 
on the left-hand side of Eq. 6.17 is called the excitation vector 


P. = - e.E. 
J J J 


(1 - a .) q». 
ej ej 


(6.18) 


on the right-hand side, the transfer matrix is given by 


where 


6 . . 
JJi 


- X- (1 - a.. ] SS. = M.. 
i 11 11 J 1 


6 . , = 
Ji 


0 i^j 

for 

1 i=j 


(6.19) 


is the Kronecker delta, ^- s a suitable quadrature coefficient and SS ^ ^ is 
given either by Eq. 6.15 or by Eq. 6.16 depending on whether j refers to the 
fin or to the channel, respectively. There is no matrix multiplication 
implied in Eq. 6.19, hence the underscores. 

In the present program phase, Eq. 6.17 is solved at every time step 
only when the transfer matrix is sufficiently temperature sensitive, otherwise 
the transfer matrix is completely inverted only once to yield the unknown 
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radiosity at any time. 


E (V 1 p j (6,20) 

i 


through a simple matrix multiplication. It may be noted that the most signifi- 
cant temperature dependence of optical properties is contained in the excita- 
tion vector P . . 

J 

4 

All radiant heat fluxes are normalized with respect to aT where 

r o 

T q designates the reference temperature, that is the fluid inlet temperature. 
Exchange factors are nondimens ional and need not be normalized. 


w = 

j 




P 



o 


( 6 . 21 ) 


The nondimens ional forms of Eqs . 6.4 and 6.20 are used to compute 

q , in Eqs. 2.17, 4.11 and 5.5. 

hiet,rad n 
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7 . Aerodynamic Heating 

The aerodynamic heating model used to evaluate the convective flux 
from the radiator surface on the orbiter vehicle is subdivided into three 
major regimes. The first regime encompasses low speed flow for which the 
heat transfer coefficients are determined from expressions appearing in stan- 
dard heat transfer texts for flow over a flat plate. The second regime con- 
sists of a model for high speed flow in which the convective heat transfer 
coefficient is evaluated from an experimental correlation for flow over the 
upper surface of the shuttle orbiter vehicle (Ref. 9) . The third section of 
the model encompasses a low to high speed transitional flow regime. Within 
this regime the heat transfer coefficient is an interpolated value that lies 
between the values obtained in the low and high speed regimes. Calculations 
for the convective heat flux for all three regimes are based on Eckert's 
reference enthalpy method (Ref. 10). 

Within each of the three regimes the heat transfer coefficient is 
calculated for cases where the flow is laminar, transitional or fully turbu- 
lent. In addition to the evaluation of the heat flux when the flow is forced, 
the procedure accounts for heat transfer by free convection at times when the 
shuttle vehicle is stationary or moving with a relatively low velocity. 

The program for the evaluation of the aerodynamic heating rate is 
divided into six sub-tasks each of which is written as a separate subprogram. 
This procedure allows for changes in the periphery of the program without 
affecting the program foundation. The basic calculations are carried out in 
and controlled from the SUBR0UTINE C0NVEC. Atmospheric temperature and speed 
of sound are evaluated within the SUBR0UTINE ATM0S. Atmospheric properties 
evaluated at the reference temperature are calculated within the SUBR0UTINE 
REFP. The orbiter velocity and altitute are evaluated in the FUNCTI0N sub- 
program ALTVEL. The SUBROUTINE NUS evaluates the Nusselt number for the 
radiator system. 

It should be noted that the analysis does not account for the effects 
of shock wave interaction or interference heating caused by flow interference 
between the orbiter, booster or any supporting structure. 
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The analysis for the determination of the aerodynamic heating first 
requires the evaluation of a reference temperature which is used for the 
determination of all air properties. The reference temperature is a function 
of the Prandtl number and recovery factor of the air, as well as the vehicle 
Mach number. 

The Mach number in turn is a function of the altitude and velocity 
of the orbiter at any instant time. Altitude and velocity profiles for the 
orbiter are contained in data arrays supplied by the user (see Par t A-3c of the User 
Manual. If the value of the Integer I is used to specify either ascent or reentry 
phase, then the N paired data points which define the velocity V and elevation 
Z as a function of time t may be expressed functionally as 

V i = v ± (I, t ± ) i = (7.1) 

Z i = Z i (I * t i ) 1 = M.-'-.N (7.2) 


Once the orbiter velocity and altitude are known as a function of time, the 
vehicle Mach number M may be calculated from the equation 

M - — (7.3) 

c 


The speed of sound c is an atmospheric property that is a function only of the 
altitude. 

The reference temperature is a function of the recovery factor r 
which for laminar flow is 



and for turbulent flow (Ref. 11) is 
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(8 + 0.528M 2 / (22 + M 2 )) 


r = N, 


Pr 


(7.4) 


where the Prandtl number N p ^ is an atmospheric property. To avoid a disconti- 
nuity in the value of the recovery factor between the laminar and turbulent 
flow models, Eq. 7.4 was used as the expression for the recovery factor for 
both flow models. The resulting error in the reference temperature was found 
to be approximately 5 R in an extreme case. 

All of the properties for the atmosphere used in the evaluation of 
the heat transfer coefficient are evaluated at the high speed reference tem- 
perature. Eckert (Ref. 10) recommends the expression 

i* = 0.5 (i„ + i w ) + 0.11 r ( Y - 1) (7.5) 


for the reference enthalpy i* which can be converted to the reference tempera- 
ture T* once the relationship 

T = T (i ) (7.6) 


between the atmospheric enthalpy and temperature is known. The subscripts 
and "w" in Eq. 7.5 refer to the enthalpy of the air evaluated at the free 
stream and surface temperature, respectively, and y is the ratio of the 
specific heats for air. 

It should be mentioned that when velocities are low (M ->• 0) Eqs . 7.5 

and 6 yield a reference temperature that approaches the film temperature 

(T + T )/2. As a result Eqs. 7.5 and 6 were used to evaluate the reference 
N 00 A 

temperature for all three flow regimes, i.e., low speed^high speed and the 
transitional regimes. 

The convective flux is the product of the convective heat transfer 
coefficient and the difference between the air enthalpy evaluated at the sur- 
face temperature and at the adiabatic wall temperature. The adiabatic wall 
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enthalpy i is related to the free stream enthalpy, recovery factor and 
vehicle velocity by the relationship 


i 

aw 


+ 



(7.7) 


The convective heat transfer coefficient h used in the reference enthalpy 
method may be expressed in terms of the Nusselt number N 

Nu 


N 


Nu 


h . xc 
- 2, P 
* 
k 


(7.8) 


where c* and k* denote the atmospheric specific heat and thermal conductivity 
evaluated at the reference temperature T*. The symbol x denotes a characteris- 
tic length of the radiator system which for forced convection is the distance 
from the stagnation point on the shuttle to the center of the radiator panel. 

The expressions for the orbiter Nusselt number selected for the low 

speed and high speed regime and for laminar, transitional and turbulent flow 

are summarized in Table 1. Values for the Nusselt number for conditions lying 

between the low and high speed regimes were obtained by interpolation so that 

the convective heat flux from the shuttle varies continually from one regime 

to the other. The symbols N and N appearing in the table denote the 

Ke JN u ' 

Reynolds number and Prandtl number, respectively, where 

P*VX 


N. 


Pr 



* 


k 


* 


The superscript on each property indicates that the property is evaluated 
at the temperature T*. 



TABLE 1 


FORCED CONVECTION NUSSELT NUMBER FOR ORBITER 
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For the low speed regime the expressions for the Nusselt number are 

* 

those for laminar transitional and turbulent flow over a flat plate (Ref. 6 ). 
Nusselt number relationships for high speed flow regime are taken from Ref. 9 
where experimental wind tunnel data are presented for a delta space shuttle 
orbiter. The data are for leeward surface heat transfer at angles of attack 
between 10° and 30° and Mach numbers of 8 and 16. The Nusselt number is 
shown to be relatively independent of angle of attack so that the high speed 
correlation may be applied to both the ascent phase for which the angle of 
attack is approximately zero and the reentry phase when the angle of attack 
may approach 60°. The scatter in the data of Ref. 9 from the selected 
Nusselt relationships is on the order of 100%. 

The leeward surface correlation was selected because the aerodynamic 
heating rates in this region of the shuttle are relatively low when compared 
to heating rates for the lower body or windward surface. Estimates place the 
peak reentry temperature of the lower surface stagnation line around 2100 F 
while the peak leeward temperature is estimated to be about 600 F (Ref. 12). 
Therefore, placing the radiator on the upper surface of the orbiter not only 
will result in a more efficient operation upon reentry, but also will minimize 
the need for reservicing of the surface coating on the radiator panels. 

The aerodynamic heating analysis includes free convection from the 
radiator surface during pre-launch operation and when the shuttle vehicle is 
moving with a relatively low velocity. The expression for the free convection 
Nusselt number is a function of the Grashof Prandtl product where the Grashof 
number N is given by 


N 


Gr 


= g3 


(% 2 

y 


(T - 
w 


T ) 


where g represents the acceleration of gravity, p,p and T are the atmo- 
spheric density, dynamic viscosity and temperature, respectively and T^ 
denotes the radiator surface temperature. The symbol y denotes the overall 
dimension of the radiator panel in the direction parallel to the acceleration 
of gravity. Since the atmosphere is assumed to be an ideal gas, the 

•k 

pp. 296 and 313 
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coefficient of thermal expansion 8 is simply the reciprocal of the average 
absolute temperature of the air 


6 = 


T 


2 

+ T 


The free convection Nusselt numbers % Uf and their applicable ranges 
of Grashof Prandtl product used for this analysis are 





= 1.585 (N N ) °.195 
Gr Pr 

, ,0.250 

= 0.590 (N Gr N pr ) 

, V 0.333 

= 0.130 (N_ N ) 

Gr Pr 


(10 _1 < N Gr N pr < 10 4 ) (7.9a) 


(10 4 <_ N Gr N pr < 10 9 ) (7.9b) 


(N Gr N Pr > 10 > 


(7.9c) 


The values for the free convection Nusselt number given in Eq. 7.9 were 
multiplied by the ratio (x/y) and then added to the forced convection Nusselt 
number to obtain a value that accounts for combined free and forced convection 
in the low speed regime. 

Equations 7.1 through 7.8 combined with the appropriate Nusselt 
number relationship from Table 1 for forced convection and Eq. 7.9 for free 
convection are sufficient to determine the convective heat flux into the 
radiator surface which is given by 


q " = h. (i - i ) 

aero l 00 aw 


(7.10) 


The convective flux may be normalized by dividing by the heat flux crx or 

o 
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aero 


h. (i - i ) 
i 00 aw 

4 

a T 


( 7 . 11 ) 


where T q is the fluid temperature at the inlet plane to the flow channel. The 
normalized convective heat flux given by Eq. 7.11 is used in the energy equa- 
tion for both the fin (Eq. 2.18) and the meteoroid protection layer (Eq . 5.5). 
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C. DESIGN PARAMETERS 


Certain design parameters are necessary for the design specifica- 
tion, for the selection of the optimum radiator system and even the system 
definition as required for the heat transfer analysis. In the following are 
discussed, in that order, the prediction of the meteoroid protection layer 
thickness, the system weight and a collection of non dimensional groups which 
define the radiator system, the operational conditions, and the performance 
characteristics . 


& 
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8. Meteoroid Protection Thickness 

In this section an engineering equation is developed to predict the 
thickness of a meteoroid protection layer required to cover all radiator sur- 
faces that might be damaged by the impact of a meteoroid. Several assumptions 
have been made during the derivation. They are: 

1. The meteoroid particle is spherical. 

2. The meteoroid flux is isotropic. 

3. Poissons distribution law describes the probability 
of an impact of a meteoroid. 

It should be mentioned that any equation used to predict meteoroid 
protection thickness is only as accurate as the experimental data used in 
that equation. Even though much information has been published in recent 
years concerning protection theories, there is still considerable question 
as to the density, velocity and mass distribution of meteoroid particles in 
outer space. In addition to these uncertainties, two basic models for 
penetration theory have been proposed within the last decade and there appears 
to be no close agreement between the two. Experimental verification of either 
model has been hampered by the fact that particle velocities used in experi- 
mental tests have only recently approached the meteoroid velocity range. In 
short, an extremely reliable theory for the prediction of protection layer 
thickness does not presently exist. 

Structural materials that can be used in this study as a protection 
layer are copper, aluminum and beryllium. While both copper and aluminum were 
selected primarily as fin and tube materials due to their superior heat trans- 
fer characteristics, beryllium was chosen for its protection capabilities. 

The penetration theory predicts a protection layer thickness that decreases 
as the modulus of elasticity increases. Therefore, beryllium becomes an 
attractive protection material due to its high modulus of elasticity and its 
relatively low density. In fact studies (Ref s . 13 and 14) have shown that 
beryllium can significantly reduce protection layer weight. 

The basic equation (Ref. 15) for the depth of penetration of a 
meteoroid particle into a target of infinite depth is 
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( 8 . 1 ) 


where 


y empirical constant generally accepted to be in the range of 
1.5 to 2.5 


V 

c 

d 

* 


density of the meteoroid particle 
density of the target material 
velocity of the meteoroid particle 
velocity of sound in the target material 
diameter of the meteoroid particle 
constant between 1/3 and 2/3 


0 constant between 1/3 and 2/3 

The ratio of the meteoroid velocity and the velocity of sound in the 
target represents a target Mach number. The velocity of sound in the target 
can be expressed in terms of the modulus of elasticity 


c = ^ 


t 


( 8 . 2 ) 


where 

E modulus of elasticity of the target material 

g^ proportionality constant relating mass units to force 
units 

P t target density 

If the meteoroid particle is assumed to be spherical, the diameter 
may be written in terms of its mass 
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(8.3) 


where 


M meteoroid mass 
P 

p meteoroid density 
P 


The probability that an exposed surface will be struck by a meteoroid 
during a period of time can only be determined after the distribution of 
meteoroids of a given mass is known. This information is usually given in 
the form of an equation such as 


F = aM p 6 (8.4) 

where 

F cumulative number of impacts of particles with mass 
M or larger per unit area per unit time 

M mass of the meteoroid particle. 

P 

The symbols a and 3 represent experimentally determined constants. Published 
(Ref. 16) values of a and 6 vary over a considerable range, but they lie 
within the limiting values , 


1.3 x 10 15 < a < 2.54 x 10 
1.11 < 3 < 1.34 

The cumulative number of impacts on a surface with a vulnerable area 
of A during a mission time of t by a meteoroid of mass of larger is then 


[ft^ day gm 
[Dimensionless] 


FAt = At a M 
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It is generally assumed (Ref. 17) that meteoroids are randomly dis- 
tributed in outer space and that each collision can be described by Poissons 
probability law. From the Poisson distribution function, the probability of 
zero events occurring P^ when the average number of events is N is 


In P = -N 
o 

Substituting the value for N gives the probability that no meteoroid of mass 
M or larger will impact on the surface of area A during time t of 


or 


In P 

o 


= At a M 

P 


M = 
P 


1/3 

aAr 
-In P 

o 


(8.5) 


To account for the fact that all meteoroids do not strike the pro- 
tection layer normally, the meteoroid velocity V may be replaced by a critical 
velocity where 


V c = V (cos X) n 


( 8 . 6 ) 


X angle between the direction of V and the normal to the protection 
surface 

n an experimentally determined constant. 

If n is selected to be unity the damage to the protection layer caused by an 
oblique collision is based on the meteoroid's normal component of velocity. 

A more conservative approach would be to set n = 0 in which case all particles 
are considered to impact normally. 
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If the meteoroid flux is assumed to be isotropic the angle depen- 
dence may be replaced by 


(cos 


A) n = 


/ 2 ..) 

\3n93+2 / 


1/33 


(8.7) 


Finally account must be taken for the fact that the meteoroid will 
not impact on an infinite target, but one with a finite thickness. As a 
result even though the meteoroid may not penetrate the protection layer, 
spalling may damage the radiator panel. To account for spalling the thick- 
ness of meteoroid protection t used should be larger than the predicted 
penetration into an infinite target or 


t = 


aP 


00 


( 8 . 8 ) 


Accepted values of a lie between 1.5 and 2.0. 

Equations 8.1 through 8.7 may now be substituted into Eq. 8.8 to 
give an expression for the meteoroid protection layer thickness. The result 
is 


t 





V 

12 (E g /p ) 
t c t 



/ — Z Yl/38 

y3n0 3 + 2 ) 


(8.9) 


where 


t - thickness of protection layer (inches) 
a - experimental constant (dimensionless) 
y - experimental constant (dimensionless) 

a - experimental constant that relates meteoroid flux to mass 
(gm/ (day ft 2 )) 

3 - experimental constant that relates meteoroid flux to mass 

(dimensionless) 
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A 

T 

P 

O 


n 


2 

~ vulnerable area requiring protection (ft ) 

- mission time (days) 

- probability of no damage caused by impact of meteoroid 
(dimensionless) 

3 

- density of meteoroid particle (gm/cm ) 

3 

- density of protection layer (lb^/ft ) 

- velocity of meteoroid (ft/sec) 

2 

- modulus of elasticity of protection material (lb^/in ) 

- 32.174 lb ft/lb r sec 2 

m r 

- experimental constant (dimensionless) 

-experimental constant (dimensionless) 

- experimental constant that describes penetration depth as 
function of angle of incident (dimensionless). 


a 


Selection of Values for Experimental Constants 


Values for the experimental constant p, 6, p and V used in Eq. 

P 

were selected from the Manned Spacecraft Center publication for meteoroid 
environment criterion (Ref. 18). 


8.9 


The values for a and 8 for meteoroids having a mass between 1 gm 

_6 

10 gm used in Eq. 8.4 are 


a = 1.888 x 10 10 gm^/ (ft^day) (8.10) 

3 = 1-213 (8.11) 


The average meteoroid density is 

3 

Pp =0.5 gm/cm 

and the average meteoroid velocity is 


( 8 . 12 ) 
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V = 20 km/sec. (8.13) 

Values chosen for the remaining constants appearing in Eq. 8.9 are 
summarized in Table 2. Values of these constants are also listed in the table 
that will yield optimistic (minimum) and pessimistic (maximum) thicknesses for 
the meteoroid protection layer. 



Recommended Value 

Pessimistic Value 

Optimistic Value 

a 

1.75 

2.0 

1.5 

Y 

1.50 

2.5 

1.5 


1/2 

1/3 

2/3 

0 

2/3 

2/3 

1/3 

n 

1.0 

0 

1.0 


TABLE 2 . Empirical Constants for Meteoroid Protection Layer Thickness 


The following is an analysis of the sensitivity that the meteoroid 
protection layer thickness has to the uncertainty in the values of the five 
parameters listed in Table 2. This information will enable the user to judge 
his selection of these constants from within the recommended ranges. 

An expression for the error in the meteoroid protection thickness 
may be obtained by taking the logarithm of Eq. 8.9 followed by differenting 
the resulting equation. This process yields the equation 


dt — da chfl 

t a y 


+ <j> 


(pAd± 

\Pt / $ 


+ 


[ 9 1 " (-=-) 


n8 , dG r n8 dn 

3n60+2 G t 3nG0+2 J n 


(8.14) 


If the symbol E^ is selected to represent 


the relative error in the meteoroid 
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protection thickness resulting from a relative uncertainty in the value for 
the parameter a, then it is evident from Eq. 8.14 that 


dt/ 1 
a da/a 


1.0 


when all other parameters are held constant. Similarly the error caused by 
an uncertainty in the value of y will be 


Ey 


dt/t 

dy/y 


1.0 


When the value for 0 is taken to be the recommended value of 2/3 
and 3 is set equal to the value fixed by MSC’s environmental model (Eq. 8.11), 
the resulting error in the meteoroid protection layer thickness due to an 
uncertainty in the value for n is 


E 

n 


dt/t 

dn/n 




nO 


3n03+2 


■) = - 0.15 


The magnitude of the error for the final two parameters <f> and 6 are 
a function of the material selected for the protection layer. In order to 
give an indication of the range of errors that can be expected for various 
protection materials, the errors were calculated for the three structural 
materials that were selected in the program: aluminum, beryllium and copper. 

If the meteoroid particle density is assumed to be fixed at the value re- 
commended by the MSC environmental model (Eq. 8.12), then an uncertainty in 
the value of <j> from the recommended value of 1/2 would cause an error in the 
protection thickness of 


= dt/t _ 

( P d<J>/ <J> 


♦in 

P-t 


which for each of the three structural materials results in the following 


errors 
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E = -0.85 aluminum 

E^ = -0.65 beryllium 

E^ = -1.44 copper . 

The high density and low modulus of elasticity of copper makes 
its protection characteristics rather undesirable. For this reason the 
error of 1.44 indicated for copper probably will never be experienced in 
practice, and this value should be considered as a limiting case. 

An uncertainty in the value of 0 from the recommended value of 
2/3 would cause an error in the protection layer thickness equal to 


E = 

e de/e 


[Bin (4 
c 


n Q 

3nQ8+2 



which for each of the three structural materials is: 

E =0.77 aluminum 

0 

E =0.20 beryllium 

0 

E =1.01 copper. 

0 

The errors calculated in the analysis are summarized in Table 3. 


TABLE 3. Ratio of Relative Error in Thickness to 

Relative Uncertainty in Empirical Constants 


Protection 

Material 

Aluminum 

Beryllium 

Copper 

Parameter 

a 

1.0 

1.0 

1.0 

7 

1.0 

1.0 

1.0 

<i> 

-0.85 

-0.65 

-1.44 

0 

0.77 

0.20 

1.01 

n 

-0.15 

-0.15 

-0.15 





61 


If the error values for copper are excluded, the protection thickness 
is most sensitive to variations in the parameters a and y and least sensitive 
to variations in the parameter n. Even though the protection thickness is 
least sensitive to the selection of n, it should be noted that the 100% 
variation between the optimistic and pessimistic value of n is the largest 
of all of the parameters. Also it should be noted that the signs on the 
values in Table 3 indicate that an increase in the parameters a, y and 6 
result in an increase in the protection layer thickness, while an increase 
in the values for <j> and n result in a decrease in the protection layer 
thickness. This fact can be verified by the choice of the values of each 
of the parameters listed in Table 2. The values labeled as those which will 
produce a pessimistic value for the protection thickness are maximum values 
for a, y and 0 and minimum values for <J> and n. 

To further evaluate the effect the meteoroid protection thickness has 
on the performance of the fin system, the temperature of the coolant fluid 
at the exit plane of the flow channel was evaluated first under the 
"pessimistic" conditions for the meteoroid protection layer, second for 
the "recommended" conditions and finally for the "optimistic" conditions. 

The results of these computer runs are shown below. 


Case 

Protection Layer 
Thickness -Inches 

Normalized Outlet 
Fluid Temperature T/T q 

Pessimistic 

0.377 

0.8855 

Recommended 

0.063 

0.8922 

Optimistic 

0.020 

0.8932 


Even though the thickness of the meteoroid protection layer varies by 
nearly a factor of 20, the resulting error in the enthalpy drop is only 

0.8932 - 0.8855 
1 - 0.8922 


. 100% = 7.15% 
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9 . The Mass of the System 

The system mass is computed, firstly as a convenience for the user 
and secondly for the purpose of the planned system optimization. The system 
mass includes 


(i) the mass of the fluid in all tubes 

d 2 n f , d 2 u 

M c = n t — J P c dz= Vc,o — 

O o 

(ii) the mass of the fins 



(9.1) 


M f = n t H L p f (s r + s t ) 


(iii) the mass of all tubes 


and 


M = s L • tt (d + s ) p 
w t w w w 

(iv) the mass of the protection layer 


(9.2) 


(9.3) 


M = n s L p [ir(d + s + s ) - & ] (9.4) 

mtmm wmr 

but it does not include the thermal coating nor the mass of the manifold 
and the fluid in the manifold. In Eqs. 9.1 through 4 represent 

n the number of tubes 

d the tube diameter 

P the density 

L the tube length 

H the fin height, distance between fin root and fin tip 

s the thickness 

while the subscripts designate p and s as follows 
c coolant fluid 


f fin 
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ra meteoroid protection layer 
r fin root 

t fin tip 

w tube wall 
o inlet condition. 


The integral in Eq. 9.1 is time- dependent and evaluated at the initial 
conditions. 
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10 . Nondimens ional System Parameters 


The governing equations in the preceding radiator system analysis 
are developed in non-dimensional form for the purpose of (i) reducing the 
number of parameters, (ii) evolving the set of relevant parameters, and 
(iii) presenting the results in a general form which is applicable to 
groups of systems rather than an individual system. A summary of parameters 
is presented here for the detailed analysis discussed in the preceeding 
chapters . 


The transient flow field is the coolant channel and the temperature 
field over the fin can be represented as functions of: 


a) the independent variables 


time 



axial coordinate 



radial coordinate 

transverse coordinate 
normal to the channel 
axis 




( 10 . 1 ) 

( 10 . 2 ) 

(10.3) 


b) 


the dependent system variables 


fin temperature 

channel temperature 

meteroid protection 
layer temperature 

coolant fluid temperature 


e f ) = f- 


% (n>c;T ) = T 


W 


0 (n , t ) = -~ 

m p 

o 


9 c (CST ) -f- 

O 


(10.4) 

(10.5) 


( 10 . 6 ) 

(10.7) 



65 


coolant fluid pressure 
coolant fluid velocity 


tf(C,T ) = 

p o 

^<C;t ) = — 

W 

o 


( 10 . 8 ) 

(10.9) 


The solution to the problem will depend on the geometry of the 

system, the material properties and the definition of operational conditions. 

There was no attempt made to establish similitude with respect to the 

material properties because the scaling laws would either be too restrictive 

to allow for general property variations or too complex (for instance, the 

concept of corresponding states for gases). Consequently, the <j> -parameters 

defined in Eqs. -3.27 through 3.30, in Eq. 4.4 and in Eqs. 5.1, 2 and 6 are 

omitted from this summary; they constitute temperature and pressure 

variation of properties. This leaves the following list of parameters, 

n 

in addition to n, the number of tubes: 
c) the geometric parameters 

fin height- to- length ratio 
fin profile slope 

fin root thickness 


H = 


H 


s - s 
r t 

2H 


r H 


( 10 . 10 ) 


( 10 . 11 ) 


( 10 . 12 ) 


tube diameter- to- length 

. J s= — 

ratio 0 L 


(10.13) 


2s 

channel wall thickness s = — ~ (10.14) 

w d 


protection layer thickness 



(10.15) 



d) the operational parameters 


dw 

coolant flow Reynolds number N - — 2 . 

Re v 

o 


(10.16) 


Prandtl number 

(representing coolant selection) 




c o 


(10.17) 


inlet pressure heat 


compressibility 


inlet coolant power flux 


F = (£-) 

P pw o 


Trnpd^c 

Q_ -a 

4oT 



where n is the number of tubes 


incident radiant heat flux 


meteoroid velocity 


Q rad 


M = 
m 


a q 
o 


e aT 4 
o o 


m 


m 


protection layer density 

(representing selection of 
protection layer material) 



(10.18) 


(10.19) * 


( 10 . 20 ) 


( 10 . 21 ) 


( 10 . 22 ) 


(10.23) 


where p is the density of the meteoroids. 
r mt 
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Similarity for ascent and reentry operations is difficult to establish 
unless one restricts oneself to similar velocity-altitude profiles which 
can be represented by the 

max. Mach number M max = ^ V//c ^max (10.24) 


and its corresponding (through same altitude) 


Reynolds number 


( V> 


p vL 


(10.25) 


Pr and 1 1 numb e r 


k 


U C 
oo poo 


(10.26) 


Grashof number 

(before launch) 


N 


Gr 


RL 3 6AT 

2 

v 


This completes the list of non-dimensional groups resulting from 
the detailed analysis. 
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D. Extension of the Analysis to Related Systems 
XI. Non-Symmetrical Heating 

The purpose of this section is to describe a segment of the program 
which is intended to extend the application of the radiator simulation to 
non- symmetrical flow conditions in adjacent tubes. This program will also 
account for, in an approximate manner, a radiator panel which does not lie 
in one plane, or a panel which consists of U-shaped tubes. 

In the original simulation of the radiator system, if: was assumed 
that the coolant fluid entering all tubes has identical properties. As a 
result, the operating conditions and heat loss is determined for a single fin 
segment between two adjacent tubes and these conditions are assumed to exist 
for all other fin segments. In reality, this situation will exist only when 
all tubes are fed by a relatively large manifold whose flow rate and tempera- 
ture entering each tube is unaffected by the removal of coolant fluid at 
each succeeding channel entrance. 

The actual flow situation could possibly be far from the idealized 
case which was assumed to exist because it leads to a simplified mathematical 
model. A manifold of realistic size will lose heat by radiation from its 
surface and by conduction into the fin elements so that the coolant fluid 
entering individual flow channels will experience a difference in temperature 
Furthermore, unless the manifold is carefully reduced in size after passing 
each tube inlet, neighboring tubes will not receive identical mass flow rates 
Also, two adjacent tubes may be fed by the outlet flow from separate fuel 
cells which may be operating under different load conditions causing non - 
symmetrical conditions in adjacent fin panels. In short, a situation where 
two tubes receive fluid at different inlet conditions is to be expected. 

For the purpose of extending the application of the main program to 
situations discussed above, a series of program units were written and inte- 
grated into the main program unit. These programs calculate the location of 
the adiabatic plane on a fin which separates two tubes having different inlet 
flow rates and temperatures. The fin height to the adiabatic plane is then 
used as input to the main program. The calculations with the main program 
remain unaffected since the mathematical model requires only that the input 
value given for the fin height is that distance from the tube to the location 
of the adiabatic plane. 
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The analysis considered here may also account for, in an approximate 
fashion, fin elements that do not lie in one plane. If the radiator panel 
is wrapped around a cylindrical structural component (Fig. 3a) adjacent fin 
panels will be under the influence of different effective sink temperatures. 

In this analysis the assumption is made that the regular breaks in the radiator 
panel occur at the adiabatic plane between the tubes rather than at the tubes 
themselves (Fig. 3b) which is probably the most reasonable location. If the 
cylindrical panel is composed of numerous tubes and the curvature is gradual, 
the difference between the two cases should be negligible. 

The program accepts two sink temperatures supplied by the user or it 
calculates the correct sink temperature from the MRI incident flux data. The 
analysis of this section assumes steady state , one-dimensional conduction in 
an untapered fin. The fin may radiate from one or both sides when the fin 
panel is in one plane. When the panel is curved the fin is assumed to radiate 
only from the convexed surface. Fluid and fin properties are assumed constant 
at the inlet conditions. No convection from the fin surface and no radiant 
blocking of the tubes is considered. The only incident fluxes considered in 
this analysis are those accounted for in the MRI program; i.e. solar, earth 
and earth albedo. 


The user is cautioned that the analysis considered here assumes a 
one dimensional model and as such represents an approximation to the actual 
operating conditions of the non- symmetrical radiator. This caution is par- 
ticularly applicable to the panel that contains U-shaped tubes, because this 
situation can lead to a highly two dimensional condition. 

The analysis considers the ith tube of a radiator panel of thickness 

t, inside tube diameter d, tube length and distance between tubes of 2H 

(Fig. 4). The energy transfered from the ith tube carrying fluid with an 

inlet temperature of T . , specific heat of c and mass flow rate of m. is 

oi p A 


q. = m. 


c 

P 


(T . - T. .) (1 
oi bi 



( 11 . 1 ) 


The quantity U is given by 


U. 

i 


TTdh . L . 
l l 


m.c 
i P 


( 11 . 2 ) 
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(a) Fin Segments Inclined At Tube Locations- 
Geometry For Physical Model. 



(b) Fin Segments Inclined At Location Of Adiabatic 
Planes-Geometry For Mathematical Model 


Fig. 3 


Location of Adiabatic Planes 



Figure 4. Fin Geometry for Non-Symmetrical Panel 
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where h is the convective heat transfer coefficient of the fluid in the tube 
1 

and Tb. is the average fluid bulk temperature (also assumed to be the fin root 
temperature. ) 


The convective heat transfer coefficient h. in Eq. 11.1 is evaluated 
from the Nusselt number which is related to the Reynold and Prandtl numbers 
by the expressions 


N. 


Nu 


5 + 0.025 (N_ N ) 
Re Pr 


0.8 


for 


N <0.1 

Pr 

N > 2300 
Re 


N = 0.023 (N 0 -Y* 3 ) for 
Nu Re Pr 


| N >0.1 

i Pr 

N„ > 2300 
Re 


N. 


Nu 


0.0668[n N_ (d/L)] 

3.65 + SL2E m N S 2300 

1 + 0. 045[N N^ (d/L )r f 
Re Pr 


(11.3 a) 
(11.3 b) 
(11.3 c) 


The heat rejected from the coolant fluid given in Eq. 11.1 under the 

assumption of steady state and one -dimensional heat transfer can be determined 

by calculating the net radiant loss from the fin surfaces. If one side of the 

ith fin surface is radiating to an environment with an effective sink tempera- 
* 

ture of T. then 

l 


4 *4 

q, = \ L,. a (Tj - T„. ) 


li 


li 


+ n 2 . e h 2 . L. a <4 - if) 


(11.4) 


where 7] and 7]^^ are the fin effectivenesses for the fin attached to the left 
and right of the ith tube, respectively, and h^ and are the distances 

between the tube centerline and the adiabatic plane foi^ the fin attached to 
the left and right of the ith tube. 

When the fin radiates from both sides, the effective sink temperature 
will be different for both sides and the net radiative flux will be 
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q. - 21L. € h,. L. a (T 4 - T* 4 ) 

M i lx li 1 bi 1 


+ 2H 2 . c h 2 . L. a (T 4 . - T* 4 ) 


(11.5) 


*4 


where the sink temperature T\ becomes 


* 

T. = 

l 


*4 *4 

T + T 
ui li 


1/4 


( 11 . 6 ) 


The subscripts u and l in Eq. 11.6 denote the sink temperatures for the upper 
and lower surfaces of the fin respectively. 

The fin effectiveness in Eqs. 11.4 and 11.5 is a function of the 
conductance parameter N c and the ratio of the sink temperature to base tem- 
perature T^/T^ or for the ith tube 


* 


- n < N cu> VV 


* 


TL . = T] (N _ . , T./T. .) 
2i 1 c2i x bi 


(11.7) 

( 11 . 8 ) 


The conductance parameter is defined as 


N 


cli 


3 2 

€ cr tf. h, . 
bi li 

Kt 


and 


(11.9) 


N 


c2i 


m 3 , 2 
€ O' T, . h_. 
bi 2i 

Kt 


( 11 . 10 ) 


where K is the thermal conductivity of the fin material. 

In addition to having a continuous slope in the fin temperature 
profile at the adiabatic plane, the fin temperature must also be continuous 
at this plane. For a one dimensional fin the temperature at the adiabatic 
plane T is a function of the conductance parameter and the ratio of sink tem- 
perature to base or 
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T = T (N , T'7T, ) (11.11) 

N c b 

As the program searches for the position of the adiabatic plane , 
the distances between the tube centerline and adiabatic plane is systematically 
varied until the fin tip temperatures given by Eq. 11.11 for two adjacent 
tubes are equal. As an example, consider the ith and i+1 tube. Using the 
subscript 1 to denote a fin segment attached to the left of a tube and the 
subscript 2 to denote the fin segment attached to the right of a tube, the 
tip temperature profile will be continuous when 


T - T — » 0 

i 2i li+1 

Xf the difference in these two tip temperatures is denoted by 6^, then 6^ 

becomes a function of T 2j , and T u+1 . By combining Eq. 11.11 and the definition 

of N these two tip temperatures may be specified by the functional relationships 
c 

T_. = f (T, h 0 .) and (11.12) 

2i bi 2i 


li+1 


= f (T 


bi+1 ’ 


h li+l^ 


(11.13) 


The fin base temperatures of the ith and i+1 tube can be determined by Eqs. 
11.7 through 11.9 or written functionally 


T bi " f (h li’ h 2i> a " d 


T bi+1 ~ £ (h li+l’ h 2i+l ) 


(11.14) 

(11.15) 


If Eqs. 11.12 through 11.15 are combined, the difference in the two 
temperatures 6. may be expressed as 


6 i “ 6 i (h li’ h 2i’ h li+l’ h 2i+l ) 


But from the geometry of the radiator system 


(11.16) 
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h 


2i-l 


and 


(11.17) 


li+1 


= 2H i + i 



(11.18) 


where 2H. and 2H. - is the distance between the centerlines of tubes i-1 and 
l l+l 

i and tubes i and i+1, respectively. 

By combining Eqs. 11.17 and 11.18 with Eq. 11.16, the final functional 
relationship for the difference in the fin temperatures at the adiabatic plane 
is given by 


6 . = 6 . 

i l 


(h 


2i-l’ 



h 2i+P 


(11.19) 


In other words 6. can be expressed entirely in terms of the distance from tube 

l 

centerline to the adiabatic plane for the tube under consideration plus the 
same distance for tubes on either side. 

The approach to the solution for the location for all adiabatic planes 
is one which gives 6 = 0 for all i tubes. The program utilizes a Newton iter- 
ation to determine an appropriate change in ith adiabatic plane Ah_^. Written 
in matrix notation, the relationship between Ah_^ and the difference in tip 
temperatures 6 is 


I« t J = iPyKAh.J 


where the elements of the matrix p. . are given by 


[P. .1 

ij 


3h. 

3 


( 11 . 20 ) 


( 11 . 21 ) 


and the elements of the vector Atu are given by 


at- u k+1 u k 
Ah. = h. - h. 

li l 


( 11 . 22 ) 
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where k denotes the k th iteration used to find a set of adiabatic plane 

locations for all fin elements which result in 6 = 0 for all i tubes. A 

mathematical subroutine is used to invert the matrix [p. . ] so that the values 

ij 

for all adiabatic planes may be calculated from 

[Ah i ] = r P;L .r 1 [6 ± ] (11.23) 

Once the distance to all adiabatic planes have been determined 

for all tubes which lead to continuous tip temperatures and continuous, 

zero slopes at the fin tips, the net heat transfer from each fin segment 

may be determined from Eqs . 11.5 through 11.10. The exit fluid temperature 

from each tube T may then be determined from 
e • 


q - m cp(T e< - T 0 .) 
x 1 . 


(11.24) 


Since the main program unit requires only a single distance 
between the tube centerline and the adiabatic plane, the two distances h^ 
and h£^ must be averaged to yield a single value. The averaging process 1 
produces a value h^, which results in the same heat transfer for the 
symmetrical fin. This single value for each tube is calculated by using 
equation (11.5) resulting in 


h. 

i 


Hi hi + n? ho 
± 1 i *1 


n l. 

l 


+ r\ 2 

i 


(11.25) 


As previously defined, the subscripts 1 and 2 denote the two 
halves of the fin attached to both sides of the ith tube. Values for the 
fin effectiveness are evaluated from Eq. 11.7. 

The program considers the case of the radiator panel with U-shaped 
tubes to be approximated by three sequential parallel tube sections. 

The program first calculates the fin base temperature and the 
distance to the adiabatic plane for the first section of the three tube 
segments. The heat transferred from the fluid is then calculated from 
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Eq. 11.1. The outlet fluid temperature from the first tube segment is 
then calculated from the expression 

1 “ ic p (T e - V 

The exit fluid temperature simply becomes the inlet fluid temperature to 
the second tube segment (the base of the U) . The same process is repeated 
for the remaining two tube segments resulting in two more values for the 
distance to the adiabatic plane. The single value for the distance to the 
adiabatic plane necessary for entry to the main program unit is determined 
from an average weighted with respect to the length of each tube segment. 
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E. Property Fundamentals 


The fundamental principles used to prepare the required thermo- 
physical properties for inclusion into the computer code are exhibited in 
this chapter while the specific details concerning the materials treated 
in this program phase are placed into the appendices. 

The principles involved are those of macroscopic thermodynamics 
treated in most elementary texts. The approach of deriving analytic ex- 
pressions for the required properties is not unique because the starting 
point is dictated by the availability of experimental data. The result, 
however, must be of the same form regardless of whether, for instance, the 
coolant fluid is gaseous or liquid. 

The properties of the structural materials are the least problematic 
ones since they depend on the temperature at most; and the standard poly- 
nomial collocation methods are entirely sufficient. Care must be exerqised, 
however, that the collocation imply continuous fourth derivatives for highest 
integration efficiency or, less desirable, at least continuous representation 
of the property itself which may exclude piecewise allocation of degrees 
higher than one. 

What is said about the properties of structural materials holds 
in principle for the description of the atmosphere whose properties depend 
only on altitude. Even though the optical properties of the thermal 
coating depend, in general, on wave-length and temperature, the spectral 
dependence is integrated into the averaged ("gray") properties (see Eqs. 

6.5 and 6), and the results are functions of one variable, the temperature. 
Consequently, there remains but the discussion of the coolant fluid proper- 
ties which depend strictly on two state variables. 

In macroscopic thermodynamics there are required two equations 
of state for the description of a substance, namely the thermal equation 
of state f(p,p,T) = 0 which relates any one of density P, pressure p and 
temperature T to the remaining two, and the caloric equation of state, 
perhaps in the form of c° = f(T) where c° is the zero-pressure specific 
heat at constant volume. These two equations are sufficient to develop all 
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of the required thermodynamic functions, namely: 

(i) specific heat at constant pressure c^(p,T) 

(ii) isobaric thermal e'kpansion coefficient B(p,T) 

(iii) isothermal bulk modulus k(p,T) 

(iv) enthalpy h(p,T) 

These functions are discussed in Section 12. 

The transport properties, namely the thermal conductivity k(p ,T) 
and the dynamic viscosity y(p,T) are correlated on the basis of residuals 
as explained in Section 13. The properties of the atmosphere are dealt 
with in Section 14 . 
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12 . Thermodynamic Properties 

The first step in developing thermodynamic properties is to secure 
a thermal equation of state 


f(p»P»T) = 0 


( 12 . 1 ) 


For almost all pure gases and air, this equation can be found in the litera- 
ture, either in the form suggested by Benedict-Webb-Rubin (virial expansion) 
or in that suggested by Beattie-Bridgeman. Both equations are explicit in p, 

P = P (P»T) (12.2) 


so that k is immediately obtained from Eq. 11.2 by implicit differentiation 


k(p,T) 



(12.3) 


which can be evaluated as <(p,T) after inversion of Eq. 11.2 into 


P = P(P,T) . (1.2.4) 

The inversion of Eq. 12.2 is facilitated by computing (8p/3p ) from 
Eq. 11.2 and subsequently applying the Newt on - Rap h son method along the 
specified isotherm with temperature T in Eq. 12.4 

The isobaric thermal expansion coefficient 



is also obtained by implicit differentiation of Eq. 12.2 while keeping the 
left-hand side constant. 

After having collocated the zero-pressure specific heat at constant 
volume; that is, c° (T) by a power polynomial in T, one obtains first the 
specific heat at constant volume 
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c y (p,T) = c ° (T) - T 



(12.5) 


and then the specific heat at constant pressure 

2 

C P (P ’ T) = °v (p ’ T) + ff- • (12.6) 

The derivative in Eq. 12.5 is obtained from Eq. 12.2; and g and < are both 
function of p or p and T. 

Finally the enthalpy h is calculated from its definition 


h (p > = ^ + u (p , T) 

where the internal energy u may be obtained by two successive integrations, 
the first one along an isotherm (over p) , and the second one along an iso- 
chore (over T) : 


T P 

u = u < p o’V + // 

T 0 


c v (p’,T , )d 'dT* + 


A- ? 


dfij (12.8) 

( P ') 2 


Liquids can be treated, in principle, as gases; except that the 
equation of state, Eq. 12.1, is rarely available. One may find, with little 
difficulty the zero-pressure isobaric expansion coefficient 3^ = c(T), and 
then represent adequately the isothermal bulk modulus by 

k(P,T) = a(T) + b (T) p • (12.9) 


Under any circumstances, one must satisfy 



( 12 . 10 ) 
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which yields, form Eq. 12.9 

B(P»T) = -a'p - y~ p 2 + c 


( 12 . 11 ) 


where primes indicate differentiation with respect to T, of the polynomials 
a(T) and b(T) in Eq. 12.9. Equations 12.9 and 11 yield for the density 


T 

[a(T)p + i b(T)p 2 - f c( T’)dT'] 

P(PT) « p(0,T o )e f . (12.12) 


The specific heats and the enthalpy are to be derived as for gases (see Eqs. 

12.6 and 7). Other possibilities are to develop k from the speed of sound 

and the ratio of specific heats; but the reader is warned not to imply 

k = 0 or c = c , unless there is sufficient evidence to support these 
P v 
assumptions. 
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13 . Transport Properties 

While the thermal conductivity k and the dynamic viscosity y of 
liquids may often times be adequately represented by functions of tempera- 
ture T alone (facilitated by polynomial collocation) , these same properties 
for gases depend on density as well. It is recognized that the difference, 
or residue 

^(P) = k(p,T) - k* (T) (13.1) 

between the thermal conductivity k(p,T) and the low-pressure thermal con- 
ductivity k*(T) depends only on the density. Similarly, for the dynamic 
viscosity 

<P 2 (p) = y(P,T) - u*(T) . (13.2) 

Hence k and V can be represented by the sum of two polynomials, one in p 
and the other in T. The residuals and ^ are published for a number of 
gases or may be developed from property data (Ref. 19 for and He). 
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14. Atmospheric Properties 

For the prediction of aerodynamic heat fluxes incident on the 
radiator system during ascent and reentry (see Section 7) the evaluation 
of the following atmospheric properties are required: 

Temperature 

Pressure 

Density 

Molecular Weight 

Speed of Sound 

Viscosity 

Thermal Conductivity 

Specific Heat at Constant Pressure 

Enthalpy 

The model for these atmospheric properties is presented in two 
sections. The first covers altitudes from sea-level to 301,000 feet. 

Within this range the molecular weight is assumed to be constant and the 
temperature variation with altitude is a sequence of connected line segments. 
The second section of the model covers altitudes above 301,000 feet where 
the molecular weight decreases linearly with altitude. For this altitude 
range the approximate polynomial expressions for density and pressure sug- 
gested in Part 4 of Ref. 20 were used. Errors between the values given by 
the approximate expression and the 1962 Standard Model are less than 5% 
over the entire range of altitudes. 

Atmospheric air is assumed to be an ideal gas for all altitudes. 
Therefore compressibility effects at low altitudes are neglected. The error 
in computed densities resulting from the ideal gas assumption may be as high 
as 0.057o for altitudes below 6 miles, but becomes less than 0.01% above 12 
miles (Ref. 20). The air is also assumed to be in hydrostatic equilibrium. 

All properties except geopotential altitude, specific heat and 
enthalpy are evaluated from expressions presented in Refs. 20 and 21. The 
expression for geopotential altitude was taken from Ref. 22, while specific 
heat and enthalpy data were taken from Ref. 23, 
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The model developed for the atmospheric properties is considered 
to be applied to altitudes tip to 100 miles and to latitudes between 30 and 
60°N. It is anticipated that atmospheric properties are not needed for 
altitudes exceeding 100 miles, because the convective heat flux from the 
fin system will be negligible at this altitude and above. 


The properties for the earth’s atmosphere are known with increased 
uncertainty as the altitude increases. In fact, the 1962 Standard Atmosphere 
(Ref. 21) consists of four regions as follows: 


0 - 20 km 

20 - 32 km 
32 - 90 km 
90 - 700 km 


Standard 

Proposed Standard 

Tentative 

Speculative 


Any uncertainty in the atmospheric properties will naturally be 
reflected as an error in the convective heat flux on the shuttle vehicle. 
Fortunately during the ascent phase of the shuttle operation the convective 
flux from the radiator system is fairly small compared to the radiative 
flux by the time the shuttle has approached altitudes for which the atmos- 
pheric properties are considered to be "speculative"; on the other hand 
during re-entry, significant convective fluxes are known to exist at alti- 
tudes above 90 km. As a result every effort should be made to revise the 
existing atmospheric property model at high altitudes as new data become 
available . 

The atmospheric model is based on several primary constants. 

The sea-level pressure, temperature, molecular weight, density and acceler- 
ation of gravity and the universal gas constant were assigned the fixed 
values of 

P = 2116.22 lbf/ft 2 
o 

T = 518.67 R 
o 

M = 28.9644 

° 3 

P = 0.07647 lbm/ft 

° 2 

g = 32.1741 ft/ sec 
o 

R* = 1545.31 ft lbf/lb mole R 
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Properties for Altitudes Less than 301,000 feet , 

a. Geopotential Altitude - H 

The state variables for air are expressed in terms of the single 
variable, the geopotential altitude 



o 


Z - 1.573126 x 10“ 7 Z 2 + 2.4656553 x 1(T 14 Z 3 
- 3.8667054 x 10“ 21 Z 4 + 6.0621354 x 1(T 28 Z 5 
- 9.5013649 x 10~ 35 Z 6 


where Z is the geometric altitude in meters, g^ is the acceleration of gravity 
at sea-level and g(Z) denotes the local acceleration of gravity. See Ref. 22 
for details . 

b. Temperature - T 

The general expression of the temperature as a function of geopoten- 
tial altitude is 

T = T b + L(H - H b ) (14.2) 


T and H, are the endpoints of straight-line segments representing T(H) and 
b b 

are listed, together with L(H) in the following table. 

c. Molecular Weight - M 

The molecular weight is constant at a value of 28.9644. 

d. Pressure - P 

Within a region where the temperature varies linearly, the ideal 
equation of state and the hydrostatic yield the following expressions for 



«b 

(km) 

L 

(K/km) 

T b 

(K) 

0 

- 6.5 

288.15 

11 

0.0 

216.65 

20 

1.0 

216.65 

32 

2.8 

228.65 

47 

0.0 

270.65 

52 

-2.0 

270.65 

61 

-4.0 

252.65 

79 

0.0 

180.65 

90 


180.65 


TABLE 4. Lapse Rate and Base Temperatures for 
Atmospheric Model 
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pressure: 



[ 


g M 
o o 

R*L 


T fc + L(H 


- V J 


(L t 0) 


(14.3) 




SqV H - H b ) ' 

R*T, 

b - 


a - q) 


(14-4) 


i.e 11 £ H £ 20 km 
47 £ H £ 52 km 
79 < H < 90 km 


The subscript "o" denotes a quantity evaluated at sea-level and the subscript 
M b u denotes a quantity evaluated at the base of one of the straight line seg- 
ments of the atmospheric model. 

c. Density - p 

The density may be calculated from the ideal equation of state once 
the temperature and pressure have been evaluated. 

MP 

P r*X (14.5) 

f . Speed of Sound - c 

The speed of sound was evaluated from the expression 

1/2 

R*T "1 

Y M J 



c 


(14.6) 



89 


For altitudes less than 301,000 feet the ratio of specific heats is 
taken to have a fixed value of 1.40. 

g. Viscosity - y 

The dynamic viscosity was evaluated from the expression 


where 


and 


y = ^ 


3/2 


T + S 


3 = 1.458 x 10 


kg 

sec m(K) 


1/2 


S = 110.4 K 


(14 . 7) 


h. Specific Heat at Constant Pressure c^ and Enthalpy i 

Values for c and i between the temperatures of 100 R and 6400 R 
p 

were taken from the standard Gas Tables (Ref. 23) and placed in the program 

in tabular form. A value of c and i at any temperature intermediate to a 

P 

pair of tabular values was determined by an interpolation routine (see 
Section HI 4). 
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Properties for Altitudes Greater than 301,000 Feet 

a. Pressure - P 

The pressure for altitudes between 301,000 and 528,000 feet is based 
on the polynomial approximation given in Part IV of Ref. 20. The pressure is 
written in terms of the sea-level pressure P q in the form 

.11 

P = P o X tA n Zn J~ 4 ! (14.8) 

' r\=o ' 

where Z is the geometric altitude and values for appear in Table 4.1 of 
Ref. 20. 

b. Density - p 

The density is written in terms of a similar polynomial 

li 

P = P o X tV n J~ 4 (14.9) 

* n=o ’ 

where values of appear in Table 4.1 of Ref. 20. 

c. Molecular Weight - M 

The molecular weight is assumed to vary linearly with altitude Z 
(see Fig. 1.2.7 in Ref. 21). The resulting expression for M is 

M = 28.9644 - 0.030949 (Z - 90) 


whre Z is in km. 

d. Temperature - T 

The temperature is calculated from the values for pressure, density 
and molecular weight indicated above from the ideal equation of state 
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T = 


PM 



9i 


e. Speed of Sound - c 

For altitudes greater than 301,000 feet the equation for the speed 
of sound is the same one as used for the lower altitudes, but the ratio of 
the specific heats is no longer assumed to be equal to 1.40. The ratio of 
the specific heats varies with the molecular weight according to the expression 


c 

Y = 

c - R /M 
P 

The remaining properties are calculated using identical expressions 
to those outlined for altitudes less than 301,000 feet. 



III. NUMERICAL TECHNIQUES 


1 • Introduction 

The analysis carried out in Chapter II lead, as far as the mathematical 
problem formulation is concerned, to three initial value problems and one 
matrix manipulation. The three initial-value problems are to establish 

(1) the initial conditions for the coolant fluid, defined by 
Eqs. 3.42 through 3.45, 

(ii) the dynamics of the coolant flow, defined by Eqs. 3.40 through 
42 and 44, 

(iii) the temperature field throughout the system, defined by Eqs. 2. 18 
through 2.22 for the fin, Eqs. 3.39 and 46 for the coolant, Eqs. 

4.6, 7 and 8 for the channel wall, Eqs. 4.7, 5.3 and 5 for the 
protection layer, or Eq. 4.11 replacing Eqs. 2.19, 4.6, 7 and 

8 and 5.3 and 5 in the case where Eq. 4.9 is satisfied. These 
equations must be supplemented by the specification of the ini- 
tial, non-dimensional temperature everywhere in the system. 

Each initial- value problem is solved by a fourth-order Runge-Kutta-Simpson 
integration discussed in Section III-2. 

The radiosity equation, Eq. 6.17 requires the matrix manipulation, namely 
either the solution of a system of linear algebraic equations, or a matrix 
inversion whenever the optical properties of the thermal control coating are 
considered temperature independent. Either task is accomplished by elementary 
row operations which transform, in a single process, the augmented coefficient 
matrix into a row-reduced echelon matrix. The reader is referred for this 
trnasformation to standard texts on linear algebra (Ref. 9). 

Additional mathematical operations are programmed as subprograms which 
may be generally applied and which are discussed in Sections III-3 through III- 7 in 
this order: an evaluation of polynomials in one variable, an Aitken inter- 

polation, first and second differentiation, definite integration and integration 
with variable upper integration limit for functions of equally spaced arguments y 
and solution to system of linear algebraic equations. 
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2. The Runge-Kutta-Simpson Integration 

Two types of initial- value problems are to be solved in this program. 
The first type includes Item (i) and (ii) mentioned in the Introduction, namely 
the fluid dynamics exclusive of the transient fluid temperature field, and in- 
volves ordinary, first-order differential equations, linear in the derivatives 
with respect to the axial distance £, that is Eqs. 3.42, 43 and 44. The equa- 
tions are solved explicitly for these derivatives so as to take on the general 
form of Eq. 2.1: 


dy i 

dT = f i < x ’ y l’ y 2’-"’ y n ) 


(2.1) 

y ± (0) = a ± , i = 1,2,..., N 

(2.2) 


Equation 2.2 constitutes the appropriate initial conditions. The other type 
of initial- values problem, mentioned as item (iii) in the Introduction, in- 
volves partial differential equations which are linear and of the first order 
in the time-derivatives; moreover, all equations, Eqs. 2.18, 3.46, 4.8 and 5.3, 
are explicit in the time derivatives. Having subdivided the radiator system 
into intervals, equally spaced in each appropriate domain (fluid, wall, fin, 
etc.), and then written the different equations corresponding to each one of 
the resulting N interior nodal points, one may discretize the spatial deriva- 
tives occurring on the right-hand sides of the partial differential equations. 
The result is a set of ordinary differential equations, with time as the inde- 
pendent variable but of a form which is identical to Eq. 2.1. Equation 2.2 
is given by the initial temperature distribution; a uniform temperature was 
chosen for the first start of the integration (a^ = a; i = 1, 2, ...,N), sub- 
sequent integrations during optimization runs are expected to start from the 
previously computed steady-state temperature distributions. The boundary con- 
ditions may be satisfied in three different ways. Either, the temperatures 
are computed directly from the finite-difference equation representing the 
boundary conditions at the end of every time step, or secondly, the boundary 
conditions may be included into the system of Eqs. 2.1 after differentiation 
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with respect to time, or lastly, an equation of the form of Eqs. 2.1 may 
be derived directly from a control volume bounded at one side by the boundary 
of interest. All three possibilities have been utilized. 

Discretization introduces obviously a truncation error; all spatial 
derivatives are represented consistently with a truncation error proportional 
to the square of the local spatial interval (see Sect. Ill 5) but higher order 
terms may be included anytime by modifying a single program unit each for the 
first and second derivatives. 

The system of Eqs. 2.1 and 2.2 is solved by a fourth-order Runge-Kutta 
integration, that is if the f^ in Eqs. 2.1 have continuous fourth-order deriva- 
tives, the time-related accuracy of the integration is of order four (Refs. 26 
and 27). Under much weaker conditions, namely uniform Lipschitz continuity. 

Ref. 26, the accuracy is still first-rorder and stability is secured. It may 
be noted that the Lipschitz continuity is also the prerequisite for unique- 
ness of the solution to Eqs. 2.1. 

An existing single-precision, floating-point Runge-Kutta- Simp son subprogram, 
SUBROUTINE RKS, written by R. Schubert at the Aerospace Corporation was used. 

Its fixed-step integration mode was employed for the integration of the fluid 
flow variables along the channel axis, while the transient temperature field 
was integrated with variable time steps, chosen automatically so as to keep 
the "truncation error" per time step below a specified limit. The absolute 
and relative errors A. and R^ are specified, by the user, for each variable 
y^, and after every Runge-Kutta integration step a Simpson integration is 
carried out over the same interval and with the intermediate derivatives as 
used in the former integration. From the difference D between the two inte- 
grations is calculated the "truncation error" measure 

D. 
x 


and if 0.75 < E then the time step DEL is divided by /l0 and the step is 
repeated, if 0.075 < E <_ 0.75 then DEL is multiplied by /l0 for the 



E = max E . = 
m l 
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subsequent step. 

All variables are set equal to their initial values in the program 
which calls RKS. During the integration RKS interacts directly with two other 
subroutines, namely DERIV and CNTRL, whose names are the first elements of 
the argument list in the call statement. The first subroutine, DERIV, serves 
to compute all N derivates dy^/dx in accordance with Eqs. 2.1. The second 
subroutine, CNTRL, controls the output during integration and the termination 
of integration. Output of current values of all variables along with impor- 
tant system parameters is provided under two different integration modes: 
general transient system simulation, MST0R = 0 in NAMELIST /RUN0PT/ , produces 
output in arbitrarily chosen, fixed time steps, DTWRTE, up to the final time 
TEND, both specified in NAMELIST /RUN0PT/ and in hours; the second mode serves 
to compute the steady-state conditions and is invoked by setting MST0R = 1 
and by specifying the number LIMWRT of time intervals DTWRTE at which output 
is desired. 

The integration under the second mode (MST0R — 1) is terminated as soon 
as the expected truncation error due to program termination is less than five 
times the specified relative error per time step, RLIMIT, that is R_^ in Eq. 
2.3. The largest truncation error associated with the j-th time step is 
anticipated on the basis of Eqs. 1.1 and 2 in Section II as follows 


6. = max 6. . = max { A^x 
3 , 


Zn 


y i>3 -.I 


, i = 1,2,..., N 


(2.4) 




The maximum is taken from all N modal points, A.. x is the current integration 
step size with index j, and y^ stands for the dy^/dx in Eqs. 2.1. 


The argument list of RKS (and RKSF) is as follows: 

(i) DERIV, name of derivative subroutine } declared as EXTERNAL 

(ii) CNTRL, name of control subroutine J in calling program 
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(iii) Y , array name*, containing the y_^ T s in Eqs . 2.2 ** 

(iv) DY, array name*, containing the dy^/dx in Eqs. 2.1 

(v) A , array name*, containing the A^’s in Eq. 2.3 ** 

(vi) R , array name*, containing the R^'s in Eq* 2.3 ** 

(vii) T , the independent variable X in Eqs. 2.1 ** 

(viii) DEL, the integration step**, DEL ^ 0 


(ix) 

N , 

(integer) the number of equations** 



(x) 

IFVD 

= 

0: 

variable step size**, see Eq. 2.3 





= 

1: 

fixed step size equal to DEL 



(xi) 

IBKP 

= 

0: 

adjust step size at most once before 

repeat , ** 



= 

1: 

adjust in accordance to Eq. 2.3 



(xii) 

NTRY 

= 

1: 

continue integration**, normal start 





o 

2: 

return from RKS 


to be changed 



= 

3: 

repeat last step with new DEL 


in 

CNTRL 



= 

4: 

restart 




(xiii) IERR = 0, normal integration 

= -1 indicates singularity when IFVD = 0 
= +1 indicates denominator vanishes in Eq. 2.3 at some 
time during integration. 

(xiv) through (xx) are array names* with which the user need not be 
concerned except YS that contains the y^s in 2.1 at 
the previous times step: DELY, PD, SD, YS, YST, DYST, YSIMP. 

The SUBR0UTINE DERIV communicates with RKS only via its argument list which 
contains, in this order, Y, DY, and T, as specified above under iii, iv, and 
vii. Here, the current values of Y and T are supplied to DERIV, and the corres- 
ponding values of DY returned by DERIV to RKS . 

The SUBR0UTINE CNTRL communicates with RKS also via its argument list. 

It contains Y, DY, DEL, T, NTRY, IFVD as specified above under iii, iv, viii, 
vii, xii and x, respectively. From the array Y are available for output all 
the results of integration. The time step may be modified to reach a specific 
time value; and by specifying NTRY one controls the integration process from 

& 

Declared in calling program as array with dimension size equal to the 
number of differential equations. 

To be specified prior to the calling statement. 
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variable to fixed step size during the integration by resetting IFVD. 

This completes the discussion of the integration of both ordinary and 
partial differential equations as they occur in the analysis developed in 
Chapter II. The discussion is deemed sufficient to enable the user to apply 
the RKS routine to other problems as well. 


« 
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3 . The Evaluation of Polynomials 
All polynomials 

, , 2 , N 

z - a + a,x +a„x + ... + a„x 
o 12 N 


(3.X) 


are carried out in a function subprogram based on the simple, efficient recur- 
rence relation 


z . , , = x • z . + a . 
l+l i l+l 


i = 0,1,..., N-l 


(3.2) 


The coefficients a, i = 0, 1, . . . ,N must be placed, in the calling program, 
into an array of dimension (N + 1) , N is an arbitrary positive integer. 

The procedure is coded as a function subprogram called P0LY(N,A,X), where 
X is the argument x in Eq. 3.1, A is the array containing M = N + 1 elements 
starting with A(l) = a^. 



99 


4. Aitken Interpolation 

Experimental data and supporting computer results which are not repre- 
sented by analytic expressions are interpolated by Aitkens interpolation 
technique (Ref. 27). An (n + l)-point Largangian interpolation is reduced 
to a sequence of 1/2 n (n + 1) linear interpolations. The interval spacing 
is arbitrary; and any number M > n of ordered pairs (x^,y^) can be supplied 
in the calling program. The n points of interpolation are spaced equally 
about the point x of interpolation. It should be noted, however that unless 
n = N or n = 2 the result y(x) is not continuous in general. Care must also 
be taken that all nodes x^,X 2 > ••• are distinct. 

The procedure is coded as a function subprogram called YINT(X,Y,M,N,P) , 
where X and Y are the names of arrays that have the same dimension M and 
contain the ordered pairs (x^,y^) , i = 1,2,..., M such that x^ < x 2 < • • • 

The number n of points used for the interpolation is specified as N, and the 
value of x at which to interpolate is supplied as P . Note that 2 <_ N <_ M 
must be satisfied. 
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5 . Numerical Differentiation 

The first and second derivative of tabulated functions of equally spaced 

arguments is carried out in SUBR0UTINE DDX(Y,DY,DX,N) and in SUBR0UTINE 

D2DX2 (Y ,D2Y ,DX,N) , respectively. Each subroutine requires that two arrays 

be declared in a DIMENSI0N statement in the calling program, fo have a£ least 

N elements ; one array for the set of ordinates Y -* y^ supplied by the calling 

program, the other array for the return of the results, that is DY -> dy./dx 
2 2 

or D2Y -> d y./dx , The argument interval Ax and the number of ordinances y. 

are to be specified as DX and N, respectively. However, in order that terms 

of order Ax be retained including at the endpoints of the domain, N must be 

no less than 3 for DDX and 4 for D2DX2. The truncation error is of order y’ 11 
2 TV 2 

(Ax) and y V (Ax) , respectively, for DDX and D2DX2 . 
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6 . 


Numerical Integration 
The definite integral 

F - J y(x i )dx 

X, 


1 < i < N; N >_ 2 


and the indefinite integral 


G. 

J 


X? 


■/ 


y(x ± ) 


dx + 


GCXf) 




1 < 1 < N ; 1 < j <N ; N > 3 


of a tabulated function y^ of an equally spaced argument, x^, 
x-^ + Ax, x^ + 2 Ax,...„,x^ + (N - l)Ax is carried out by a modified Simpson 
integration in the FUNCTI0N DEFINT (Y,DX,N) Subprogram and in the SUBROUTINE 
FINT(Y,YO,DX,N,F) , respectively. 


For DEFINT the ordinates y^ are to be placed in the array Y whose dimen- 
sion of no less than N elements must be declared in the calling program. The 
argument interval and the number of ordinates are specified as DX and N, 
respectively. 

For FINT there are two array declarations necessary in the calling pro- 
gram, both for at least N elements; one for the integrand Y -> y^ and the 
other for the integral F -> G^. The integration constant G(x^) , the argument 
interval and the number of ordinates are to be supplied as Y0, DX and N, 
respectively. 

The truncation error of composite Simpson integration is 


1 18q Xn (Ax ) 4 y (lv) (5) with X]L < (?) < 
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7 . Solution to Systems of Linear Algebraic Equations 

Systems of linear algebraic equations are solved by the Gauss - 
Jordan elimination process. The same technique is used to invert matrices. 

Consider the system of n linear algebraic equations 


n 


j-1 


A. .X. = y. 
i] ] i 


The solution is obtained by performing that sequence of elementary 
row operations on the augmented coefficient matrix 



which leads to the row-reduced echelon matrix 



Elementary row operations are defined by 

(i) multiplication of a row by the scalar c ^ 0 

(ii) replacement of the r-th row by the r-th row plus c times 
the s-th row; c^O, r^s; r,s <. n 

(iii) interchange of any two rows. 
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The augmented coefficient matrix consists of the coefficient matrix A in the 

ij 

first n columns and the known column vector y^ in its last, (n+l)-st column. 
The row-reduced echelon matrix has the properties that 


(i) the first non-zero element in each non-zero row is 1, 

(ii) every zero-row occurs below every non-zero row, 

(iii) if the first r rows, i = 1,2, ..., r have their non-zero 

entry in column k. then the k. f s satisfy k- < k„ < . . . < k . 

I l 12 r 

In our specific case, the row- reduced echelon matrix has the identity matrix 
in place of the coefficient matrix. 

When a matrix is inverted then the augmented coefficient matrix 
consists of the n x n coefficient matrix in its first n columns and the n x n 
identity matrix in the second n columns, j = n+1, n+2, ..., 2n. The process 
indicated above leads n x n identity matrix in the first n columns and 
the inverted coefficient matrix in the second n columns, from j = n+1 through 
j = 2n. 


The particular elementary row operations required are 

(i) division of the i-th row by A^, 

(ii) subsequent multiplication of the resulting i-th row by the 
element of the k-th row, 

(iii) and subsequent replacement of the k-th row by the difference 
between the k-th row and the i-th row. 


This process has to be repeated, essentially, for each row i = 1, ..., 


n. 
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IV RECOMMENDATIONS AND CONCLUSIONS 

The analysis described in this report has been successful in 
simulating the transient heat transfer characteristics of a radiator 
system under operational conditions expected in flight. The analysis 
serves as a basis for a rigorous computer program that has been system- 
atically sub-divided into modular subprograms. The modular concept 
facilitates repeated simulation runs with different structural materials, 
meteoroid protection material, coolant fluids, and thermal control 
coatings . 

The program predicts the net system heat rejection, the fin, tube 
and fluid temperature profiles, fluid pressure field as well as the 
meteoroid protection layer thickness and mass of the entire system. 
Optimization of the system performance may be achieved through enumeration 
of the parameter sets. 

The program has been thoroughly checked and run for a large number 
of different simulation cases. A sample representation of these cases 
may be found in reference [29]. These sample runs have aided in the 
design of the radiator system during ascent, reentry, transient orbital 
and steady state orbital conditions. As a result of these computer 
runs it is recommended that an unprotected radiator not be used during 
reentry phases of the shuttle operation, because during reentry the 
aerodynamic heating has been found to exceed the ability of the radiator 
to reject heat. Also experience gained from the sample runs has shown 
that the rigorous analysis should not be used for parameter studies 
of the heat rejection system until a satisfactory optimum domain has 
been identified by the use of the Simplified Analysis 130]. 

Displays of typical results are presented in that report. 



APPENDIX A 


STRUCTURAL MATERIAL PROPERTIES 
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The thermodynamic and mechanical properties of all materials which 
make up the radiator system are summarized in Appendices A through C. 

Appendix A contains properties for the three structural materials: 
copper, aluminum and beryllium. Copper and aluminum are intended to be used 
primarily as fin and tube material, while beryllium was selected for a 
meteoroid protection material. Four properties are evaluated for each 
material: specific heat at constant pressure, thermal conductivity, modulus 

of elasticity, while (1/k) (dk/dT) is computed by differentiating the thermal 
conductivity relationship with respect to temperature. 

Appendix B contains properties for four coolant fluids: helium, 

Dow Corning 200 Silicon oil, the liquid metal NaK, and two 3-M Company 
f luorochemical liquids FC-43 and FC-75. Six properties are evaluated for 
each coolant fluid: isobaric thermal expansion coefficient, isothermal 

compressibility, specific heat at constant pressure, enthalpy, thermal 
conductivity and dynamic viscosity. In addition, two equations of state are 
included for each fluid, one explicitly in density and one explicit in 
pressure, while the property (1/k) (dk/dT) is computed by differentiating 

the thermal conductivity relationship with respect to temperature. 

* 

Appendix C contains the total hemispherical emittance and two 
auxiliary radiative properties used in the program for the surface coating 
Z-93. 

All property relationships listed in these Appendices are presented 
in analytical form obtained by fitting a power polynomial through the data 
points. The data points listed in the tables are taken from the reference 
entered before each table. Numerical techniques used for the curve fitting 
process are explained in Section III. 

The polynomial expression for each property has been compared with 
the referenced data and within the listed temperature range has been found 
to deviate by no more than the percentage error indicated. 
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I . COPPER 


1. Specific Heat 

Reference: Touloukian, Y. S., "Thermophysical Properties of High 

Temperature Solid Materials," Thermophysical Properties Research 
Center, Purdue University, Vol. 1, 1967, pp . 456-7. 

Data Points : 


T 

C 

P 

600 R 

0.0920 Btu/ (lbm R) 

1000 

0.0975 

2000 

0.1112 


Polynomial Fit: 

Temperature Range : 400 to 2000 R 

Equation: 

C = (0.08375 + 1.375 x 10“ 2 * * 5 TR _1 ) x 32.174 Btu/(slug R) (A.l) 
P 

Maximum Error: There was no difference between the computed value 

and the input data within the accuracy of the computer. 

2. Thermal Conductivity 

Reference: Touloukian, Y. S., "Thermophysical Properties of High 

Temperature Solid Materials," Thermo physical Properties Research 

Center, Purdue University, Vol. 1, 1967, pp. 458-9. 
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Data Points: 


T 


k 


600 R 
800 
1000 
1200 
1400 


228.369 Btu/ (hr ft R) 

225.708 

222.805 

219.418 

215.306 


Polynomial Fit: 

Temperature Range: 500 to 1800 R 

Equation: 

k = (228.369 - 2.62067 6 - 0.04033 6 3 ) Btu/ (hr ft R) 


where 


T - 600. OR 
200. OR 


Maximum Error : 0 . 87% 


(A. 2) 


(A. 3) 


3. 


Temperature Variation of Thermal Conductivity 

Eq. A. 2 was differentiated with respect to temperature to yield 


i dk 

k dT 


(-2.62067 - 0.121 6 ) 


200 3 

(228.369 - 2.62067 6 - 0.04033 6°) 


R 


-1 


(A. 4) 


4. Modulus Elasticity: 

Reference: "Material Manual," TRW Equipment Laboratories, February 

1966, Report No. ER-6756, Contract No. NAS 9-4884, Fig. 50. 




Data Points: 



Polynomial Fit: 

Temperature Range : 500 to 1600 R 

Equation: 

Y = (16.55 - 0.4933 0 - 1.935 0 2 + 0.2283 0 3 ) 

X 1.44 X 10 8 (lb f/ft 2 ) 

where 


T - 459.67 R 
200 R 


Maximum Error: 0.447 o 



no 


II. ALUMINUM 7075 


1. Specific Heat 

Reference: Touloukian, Y. S., "Thermophysical Properties of High 

Temperature Solid Materials," Thermophysical Properties 
Research Center, Purdue University, Vol. 2-11, 1967, pp. 810-11. 

Data Points: 


T 

c 

P 

400 R 

0.182 Btu/ (lbm R) 

600 

0.209 

800 

0.226 

1000 

0.244 

1200 

0.270 


Polynomial Fit: 

Temperature Range: 300 to 1200 R 

Equation: 


c = (0.182 + 0.03616 0 - 0.011417 9 2 + 0.00233 0 3 - 0.000083 0 4 ) 

P 

(A. 7) 

X 32.174 Btu/ (slug R) 


where 


T - 400 R 
200 R 


(A. 8) 


Maximum Error : 0 . 34% 
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2. Thermal Conductivity 

Reference: Touloukian, Y. S., "Thermo physical Properties of High 

Temperature Solid Materials," Thermophysical Properties Research 
Center, Purdue University, Vol. 2-11, pp. 812-13. 

Data Points: 


T 

k 

400 R 

88.50 Btu/ (hr ft R) 

600 

100.395 

800 

105.96 

1000 

104.024 

1200 

99.18 


Polynomial Fit: 

Temperature Range: 300 to 1200 R 

Equation: 

k = (88.5 + 13.0665 0 +0.33275 0 2 - 1.758 0 3 + 0.25375 0 4 ) 

Btu/ (hr ft R) (A. 9) 


where 


T - 400 R 
200 R 


(A. 10) 


Maximum Error: 0.967 o 
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3. Temperature Variation, of Thermal Conductivity 

Equation A. 9 was differentiated with respect to temperature to 
yield 


1_ dk_ 

k dT (A. 11) 


(13.0665 + 0,6655 0 


5.25 0 2 + 1.015 e 3 ) 


200 (88.5 + 13.0665 6 + 0.33275 6 2 - 1.758 0 3 + 0.25375 0 4 ) 


R 


-1 


4. Modulus of Elasticity 

Reference: "Material Manual,” TRW Equipment Laboratories, February 

1966, Report ER-6756, NAS 9-4884, Fig. 50. 

Data Points: 


t 

Y 

0 F 

10.71 x 10 6 lbf/in 2 

200 

9.90 

400 

8.50 

600 

6.15 


Polynomial Fit: 

Temperature Range: 500 to 1200 R 

Equation : 

Y = (10.71 - 0.63 0 - 0.115 0 2 - 0.06 0 3 ) 
X 1.44 X 10 8 lbf/ft 2 


(A. 12) 



where 


T - 459.67 R 
200 R 


(A. 13) 


Maximum Error: 0.38% 



III. BERYLLIUM (1/2 - 3% Be 0) 


Specific Heat 

Reference: Touloukian, Y. S., "Thermophysical Properties of High 

Temperature Solid Materials," Thermophysical Properties Research 
Center, Purdue University, Vol. 6-11, 1967, pp . 753-4. 

Data Points: 


T 

c 

P 

800 R 

0.536 Btu/ (Ibm R) 

1000 

0.585 

1200 

0.622 

1400 

0.652 

1600 

0.680 


Polynomial Fit: 

Temperature Range: 400 to 1700 R 


Equation : 


c = (0.536 + 0.05667 6 - 0.0085 6 2 + 0.00083 6 3 ) 
P 


(A. 14) 


X 32.174 Btu/ (slug R) 


where 


T - 800 R 
200 R 


(A. 15) 


Maximum Error: 0.88% 
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2. Thermal Conductivity 

Reference: Touloukian, Y. S., "Thermophysical Properties of High 

Temperature Solid Materials," Thermo physical Properties Research 
Center, Purdue University, Vol. 6-11, 1967, pp. 757-9. 

Data Points: 


T 

k 

400 R 

108.863 Btu/ (hr ft R) 

600 

98.944 

800 

89.751 

1000 

80.80 

1200 

72.091 


Polynomial Fit: 

Temperature Range: 400 to 1700 R 

Equation: 


k = (108.863 - 10.5643 9 + 0.82683 6 2 - 0.20167 6 3 
+ 0.020167 9 4 ) Btu/ (hr ft R) 


(A. 16) 


where 


T - 400 R 
200 R 


(A. 17) 


Maximum Error : 0 . 907, 

3. Temperature Variation of Thermal Conductivity 

Equation A. 15 was differentiated with respect to temperature to 
yield 
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1 dk _ 
k dT 

(A. 18) 


1 

200 


(-10.5643 + 1.65367 6 - 0.60501 6 2 + 0.080668 6 3 ) 
(108.863 - 10.5643 6 + 0.826830 6 2 - 0.20167 6 3 + 0.020167 0 4 ) 



4, Modulus of Elasticity 

Reference: "Material Manual," TRW Equipment Laboratories, February 

1966, Report ER-6756, Contract No „ NAS 9-4884, Fig. 51. 

Data Points: 


t 

Y 

0 F 

44.36 x 10 6 lbf/in 2 

400 

40.41 

800 

33.95 

1200 

21.80 


Polynomial Fit: 

Temperature Range: 500 to 1700 R 

Equation: 

Y = (44.36 - 3.755 9 + 0.335 0 2 - 0.53 0 3 ) 
X 1.44 X 10 8 lbf/ft 2 

where 


T - 459.67 R 
400 R 


(A. 19) 


(A. 20) 


Maximum Error: 0.287. 




APPENDIX B 


COOLANT FLUID PROPERTIES 
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I. HELIUM 


1. Equation of State Explicit in Pressure 

Reference: Akin, S. W., Trans. ASME, Vol. 72, p. 751, 1950. 

This reference was used for all Helium properties and for brevity it 
is not repeated as the reference for the properties listed below. 

Equation: The National Bureau of Standards has published a Bendict- 

Webb-Rubin equation for helium; this equation was found to be 

. 2 

valid only up to the specified pressure limit of 3000 lbf/in . 
Preference was therefore given to the following Beattie -Bridgeman 
equation : 

p = P 2 [RT(1 - a)(i + B 1 ) - A] (B.X) 

, C 

where a = — x P 

T 

A = A 1 (l - a p) 

The values of the constants in Eq. B.l, in MKSA units, are: 

R = 2.07702 x 10 3 Nm/kgk 

2 4,2 

A 1 = 1.369595 x 10 Nm /kg 
B x = 3.5002295 x 10 _3 m 3 /kg 
C = 1.0000658 x 10 3 km 3 /kg 
a = 1.496103 x 10“ 2 m 3 /kg 

Temperature Range: 160 to 860 R 

Pressure Range: 2116 to 360000 lbf/f t 2 

Maximum Error: 0.095% 

2. Equation of State Explicit in Density 

Since the equation of state is needed explicit in density, Eq. B.l 
was solved using Newton-Raphson iteration method along an isotherm to 
give 
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p - pCp^ 

P i+1 = P i ' (3p/9p) T 


(B.2) 


Using Eq. B.l, one obtains 



RT + 2 (RB 1 T - A 1 



+ 3(A x a - 


CRB 


lx 2 
-)p 


(B.3) 


3» Isobaric Thermal Expansion Coefficient 

The isobaric thermal expansion coefficient is defined by the 
equation 

6 ’ - F (B ‘ 4) 

Since the equation of state (Eq. B.l) is explicit in the pressure, 
one can write (3 as : 


or 


, Op/3T) 

± P 

p Op/8p) T 


(B.5) 


3 = 


An A 2 CSm /x 
R[p + (B + Ar )p + — ^ P J ] 

T T J 


CR 2 CRB.. 

p[RT + 2p(RB 1 T - A 1 - ~ ) + 3p Z (A x a — )] 


(B.6) 


4. Isothermal Compressibility 


The isothermal compressibility is defined by the equation 
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K = I f i£.' i 

P 


(B.7) 


Making use of Eq. B.3, the isothermal compressibility can be written as 


K = l/p[RT + 2p(RB 1 T ~ A 1 - ~ ) + 3p 2 (A 1 a - ] (B,8) 


CRB, 


5. Specific Heat at Constant Pressure : 

Equation: Experimental and quantum statisticq.1 data for helium show 

that at zero-pressure, the specific heat at constant volume is 
independent of temperature 



(B .9) 


Using Maxwell's equations, the following expression is obtained 



(B. 10) 


From the equation of state, (Eq. B.l) the integration is carried 
out in closed form to give 


c v = R[ | + 6a (1 + | Bp ] 


(B . 11) 


The relation between specific heat at constant pressure and that 
at constant volume is given by: 


c = 
P 


c + 
v 


2 

TB 

pK 


(B.12) 
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Temperature Range : 180 to 900 R 

Pressure Range: 2116 to 216000 lbf /f t 

Maximum Error: 0.587 o for cp 

6. Enthalpy 

The variation of internal energy with both temperature and density 
is 

du = c v dT + [p - T (|E) ] . 

P P 


Substitution of c v from Eq. B.ll and equation of state data from Eq. 
B.l followed by integration along an isochore and an isotherm, one 
obtains 


g 

u = u + R [| (T - I ) + 3 P C (1 + + + V 

o2 o 2 t 2 t 2 t 2 1 


1 ^ ^ ^ ®i 2 2 

(p o “ p) + 2 ( — y 1 “ A i a)(p 0 “ p } 


where u q - 3.992 x 10 4 j/kg 

T = 10.938889 K 

0 3 

p Q = 4.669193 kg/in 

with T in K and u in j/kg. 


The enthalpy may then be determined from the equation 



7. Thermal Conductivity 


Data Points: 
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Polynomial Fit: 

Temperature Range: 160 to 860 R 

Equation: 


k = (0.0404 + 0.0302 6 - 0.0033 0 2 + 0.0003 6 3 ) Btu/(hr ft R) (B.13) 


where 


T - 160 R 
200 R 


Maximum Error: 0.547 o 

8. Temperature Variation of Thermal Conductivity 

Eq. B. 13 was differentiated with respect to temperature to yield 


1 dk 
k dT 


(0.0302 - 0.0066 6 + 0.0009 9 ) 


900 2 2 r 

(0.0404 + 0.0302 6 - 0.0033 6 + 0.0003 9 ) 


(B. 14) 


9. Dynamic Viscosity 

Equation: Viscosity correlations are usually based on the concept of 

residual viscosity: 

^(p) = y(p,T) - y*(T) 


(B. 15) 
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where 

is the residual viscosity (function of density alone) . 

£ 

y is the dynamic viscosity at atmospheric pressure. 

For helium, the dynamic viscosity is given by: 


y = y* = (2.58394 x 10" 5 T/°R)°* 647 slug/(ft hr) 
Temperature Range: 160 to 660 R 


(B. 16) 


Maximum Error: 0.29%. 
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IX. SILICON OIL 


The following properties are for Dow Corning 200 Silicon Oil (1 Centistoke 
at 77 F) . 


1. Isothermal Compressibility 

Reference: Gunst, S. B., "Density-Pressure Relationships for Two 

Low- Viscosity Dimethyl Siloxanes," Trans. ASME 72 , May 1950, 
pp. 401-7. 

Data Points: Variation of < with temperature at 0 psig and 500 

psig are given below: 


500 


100 F 

r 0 

12.35 x 10 in /Ibf 

11.60 x 10 

150 

16.05 

14.94 

200 

20.45 

18.82 

250 

26.25 

23.86 

300 

36.55 

31.88 


" 6 in 2 /lbf 


Polynomial Fit: 

Temperature Range: 560 to 760 R 

2 

Pressure Range: 2116 to 74116 lbf/ft 

Equation: The variation of k q with temperature at 0 psig is 

given by 


k = (12.35 + 2.9833 0 + 1.1 0 2 - 0.48333 6 3 + 0.1 0 4 )x 10~ 6 
o 


(B. 17) 


in /Ibf 


where 


0 = 


T - 559.67 R 
50 R 


(B .17a) 


k is assumed to vary linearly on the above range of pressure, 
hence 
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k = a + b p (B. 18) 

where 

a = k 

o 

(B.18a) 

_9k _ *500 - K o 

3p T “ 500 psi 

2 

and p is the pressure in psig and k is in in /lbf 

Fitting a power polynomial through ( with the same values 

for temperature as indicated in the table results in the equation 

b = (-1.5 - 0.0133 0 - 1.18 0 2 + 0.57333 0 3 - 0.1 0 4 ) x 10~ 9 

(B.19) 

in 4 /lbf 2 


where 0 is given in Eq. B.17a. 

Maximum Error: There was no difference between the computed 

and the input data within the accuracy of the computation. 

2, Equation of State Explicit in Density 

Reference: Gunst, S. B., "Density-Pressure Relationships for 

Two Low- Viscosity Dimethyl Siloxanes, "Trans. ASME, May 1950, 
pp. 401-7. 

Data Points: Values for the variation of density with temperature 

at 0 psig are given below: 


t 

p o 

150 F 

0.7767 gm/cm 3 

200 

0.7479 

250 

0.7188 

300 

0.6900 
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Polynomial Fit: 

Temperature Range: 540 to 760 R 

Pressure Range: 2116 to 146116 Ibf /f t 

Equation: 


where 


p = (0.7767 - 0.0288 6) x 1.94 slug/ft' 


T - 609.67 R 
50 R 


(B. 20) 


(B . 20a) 


The variation of density with pressure is given by 


a z - ^ = 


- 6 dT + < dp 


(B . 21) 


Integration along the isotherm T^ = 609.67 R and from 
p' = 0 psig to p’ = p, using Eq . B.18 results in 


z(p,T ) = a(T )p + b(T ) 

’ o o o 2 


Integration along the isobar p, from T' = T , to T' - T, 
using Eq. B.25, yields 

T £ T T 

z(p,T)-z(p,T o ) = p J a’(T f )dT' + b’(T')dT T - J c(T’)dT’ 

T o T o 

which after simplification reduces to 

2 


p = Pq e 


ap + ^ b p" 


(B . 22) 


where values for a and b are given in Eq. B.18a. 
Maximum Error: 0.1287o 


3. Equation of State Explicit in Pressure 

Since the equation of state is needed explicit in pressure, 
Eq. B .22 was rearranged to yield 
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a +Va + 2b In 


(B. 23) 


4. Isobaric Thermal Expansion Coefficient 

The zero-pressure isobaric thermal expansion coefficient can be 
written as: 



where 


is given by Eq. B.20. 


or 

0.0288 1 

o 50(0.7767 - 0.0288 0) R 


(B. 24) 


where 9 is given by Eq. B.20 a. 

Making use of Eq. B.21, together with the principle of an exact 
differential, one may write 



From Eq. B.18, the variation of {3 with both pressure and 
temperature is given by 

2 

$ = 8 q - a' p - y— p 2 (B. 25) 


where the prime superscript indicates differentiation with respect 
to temperature. The expressions for a and b as a function of 
temperature are given in Eqs . B.17 and 19, respectively. 
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5. Specific Heat at Constant Pressure 

Reference: Dow Corning, Bulletin 05-145, February 1966. 

Data Points: The available data for the variation of zero pressure 

specific heat at constant pressure for 2 centistokes are: 


t 

0 

c 

P 


80 F 

0.448 

Btu/ (lbm F) 

160 

0,454 


240 

0.463 


320 

0.476 


400 

0.491 



Polynomial Fit: 

Temperature Range: 540 to 860 R 

Pressure Range: 2116 to 74116 lbf/ft 2 

Equation: The above data for 2 centistokes silicon oil were 

multiplied by the ratio of c° for 1 centistoke to 2 centistokes 
at 77°F, to give the following expression for the zero-pressure 
specific heat at constant pressure for 1 centistoke silicon oil 


c° = (0.46 + 0.00471 9 + 0.00141 0 2 + 0.000043 0 3 ) 
P 

x 32.174 Btu/ (slug R) 

where 


T - 539.67 R 
80 R 


The variation of c 


P 


with pressure is given by 


c 

P 




dp' 


(B . 26) 


(B. 26a) 


(B . 27) 


Since the exponent in Eq. B.22 is small, the equation for the density, 
when expanded in a power series, may be truncated after the second 




129 


term in the expansion. Applying Eq. B.27 to the two- term 
expression for density as a function of pressure and 
temperature by integrating along an isotherm T from p' = 0 psig 
to p' = p, results in the expression for 


c 

P 



where 


(B. 28) 


I = 


Z-jP + 


ZoP + Z oP + Z,p + Zcp' 


and 


P' 2 

•i = 2(~) 

1 K 0 } 


z 2 = i t(2 r a ’ - a,,) - a2 i ] 

o 

l r T 2 b" , p o , /0 P o , „x ,1/2 , v 

z 3 = b’) - a(2 a - a") + jCa - b)z 1 3 

p o P o 

z, » ha’b' - a(a' 2 - ^ a’) + j(a 2 - b)(2 ^ a' - a") 

4 4 2 p 2 p 

o o 

+ \ abz^ 

1 1_ t 2 ry ry V ft P 

z 5 = lotV - 2a a ' b ' + (a - b) (a - T + — 1 1 


+ ab(2 — a' - a") + \ b 2 z, ] 

p z ± 

O 


where the prime superscript indicates differentiation with respect 
to temperature. The symbol a represents the zero pressure isothermal 
compressibility defined in equation B.17 and b is defined by equation 
B. 19 . 

Maximum Error: For zero pressure specific heat at constant 

pressure the maximum error was 0.065%. For higher pressures no 
experimental data were available for comparison. However, the 



130 


expression for enthalpy was numerically differentiated with 
respect to temperature at constant pressure and compared with 
the computed values of specific heat at constant pressure „ 

The comparison showed no difference within the accuracy of 
the computation. 

6. Enthalpy 

The variation of enthalpy with both pressure and temperature is 
given by: 

dh = c° dT + ^ [1 - T9] dp 

This expression was integrated along the isobar p = 0 psig from 
T' - 339.65 R to T 1 = T, and then along an isotherm T from p' = 0 psig 
to p ^p, to give 

h = { 80(0.46 + 0.00471 9 + 0.00141 0 2 + 0.000043 9 3 ) 

+ P .337.37 [ " ^ bb ' T P 5 - | (a b’ + a' b) Tp 4 
o 

1 p ' o 1 

+ ± (T (b 1 - a a’) - b(l + T f - ) )p J + - (a ' T - a (1 

o O 

p 1 9 p 1 

+ T — ) ) P + (1 + T - 2 -) p] X 32.174 } Btu/ slug (B.29) 
P o O 

where 9 is given by Eq. B.26a 
T is temperature in R 
a is given by Eq. B.17 in in 2 /lbf 
is given by Eq. B.20 in slug/ft 3 
b is given by Eq, B.19 in in 4 /lbf 2 
p is pressure in lbf/in 2 

and primes denote differentiation with temperature. 
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7. Dynamic Viscosity 

Reference: Dow Corning, Bulletin 05-153, July 1966. 

Data Points : 


t 


V 

0 F 

1.98 

centistokes 

100 

0.874 


200 

0.56 


300 

0.41 



Polynomial Fit: 

Temperature range: 460 to 760 R 

Equation: 

In v - (0.683 - 1.0845 0 + 0.3065 0 2 - 0.04 0 3 ) (B.30) 


where 0 = 


T - 459.67 R 
100.0 R 


where v is in centistokes and the dynamic viscosity is 
given by 

y = v p (b. 31) 

Maximum Error: 0.77% 


8. Thermal Conductivity 

Reference: Dow Corning, Bulletin 05-145, February 1966. 

Data Points: The available data for the variation of thermal con 

ductivity with temperature for 2 centistokes are: 




1 32 


1 t 

L- . - 

k 


1 — 

i 

! - 100 F 

0.0674 

Btu/ (hr ft F) 

j 100 

0.0626 


s 

I 300 

0.0578 



Polynomial Fit: 

Temperature Range: 360 to 860 R 

Equation: The procedure that was used for specific heat at constant 

pressure was followed to get an expression for the 
variation of thermal conductivity with temperature for 
1 centistoke silicon oil. The resulting expression for 
the thermal conductivity is given by 

k= (7.0052 - 2.2105 X 10 _3 T/R) x 10" 2 Btu/ (hr f t R) (B.32) 

Maximum Error: There was no difference between the computed 

and the input data within the. accuracy of the computation. 

9. Temperature Variation of Thermal Conductivity 

Eq. B.32 was differentiated with respect to temperature to 

yield 


I dk 

k dT 


- 2.2105x10 

7.0052 - 2.2105x10" 3 T/°R) 


i/r 


(B. 33) 
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III NAK - (78.6% K) 

The following physical properties of NaK(78.6 wt% K) were extracted 
from the latest version of the "Liquid Metals Handbook, Sodium and NaK 
Supplement" (to be published). Some typical properties of NaK are: 

Melting Point: 92 F 

Boiling Point: 1445 F 

Surface Tension: 0.00739 lbf/ft at Melting point 

Since all property values were taken from this single reference, 

the reference is omitted in each section below. 

1. Equation of State Explicit in Density 

Data Points: Values for the variation of density with temperature at 

zero pressure are given by 


t 

p o 

200 F 

53.21 lbm/ft 3 

500 

50.68 

800 

48.15 

1100 

45.62 

1400 

43.09 


Polynomial Fit: 

Temperature Range: 660 to 1860 R 

Equation: 

P = (58.773064 - 0.008433 T/R)/32.174 slug/ft 3 (B.34) 

o 

Since the isothermal compressibility for NaK is assumed to be 
independent of pressure, Eq. B22 reduces to 


P 



e 


Kp 


(B.35) 
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where p is the gage pressure. Since the exponent is small, the 
power series expansion for Eq . B.35 may be truncated after the 
second term and the variation of density with both pressure and 
temperature is given by 

P = (1 + k p ) (B.36) 


An expression for k is given later in this section. 

Maximum Error: At zero pressure, there was no difference between the 

computed and the input data, within the accuracy of the computation. 
For higher pressures there were no experimental data available for 
comparison. 


2. Equation of State Explicit in Pressure 

Since the equation of state is needed explicit in pressure, Eq. B.36 
was rearranged to yield 


1 P 

P * 7 (o - - D 


Isobaric Thermal Expansion Coefficient 


The isobaric thermal expansion coefficient is defined in Eq . B.4. 
From Eq. B.34, the zero pressure isobaric thermal expansion coefficient 
is given by : 


6 = 
o 


0 o 008433 I 

(58.773064 - 0 .00843 3TR“ -*-) R 


(B . 38) 


Due to the lack of experimental data, the isobaric thermal 
expansion coefficient was assumed to be independent of pressure. 
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4. Isothermal Compressibility 

In view of the experimental difficulties associated with the 
measurement of isothermal compressibility at elevated temperatures, 
such data are not generally available for liquid metals. However, the 
well-known relationship between velocity of sound c, density P, and 
isentropic compressibility is 

l_ 

K s p c 2 (B .39) 


which makes an alternative approach to the problem possible, if 
velocities of sound can be measured. Under these circumstances the 
isothermal compressibility may be obtained from the relation 


K = YK 
S 


(B.40) 


where 


c 



The relation between c and c is given by 

p v 


c 

P 



From Eqs. B.40 and B.41, one gets 


K 


K 

S 



or 


(B.41) 




K 

O 


(B.42) 
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Due to the lack of experimental data Eq. B.42 was evaluated at the 

absolute pressure p Q of one atmosphere and the temperature of 

T = 1260 R. 
o 

at T = 1260 R 
o 

P = 48.15 lbm/ft 3 
o 

c = 7544 ft/sec 
o 

6 = 1.75149 x 10" 4 1/R 

o 

C = 0.2091 Btu/ (lbm R) 
p,o 

5. Specific Heat at Constant Pressure 

Data Points: For zero-pressure specific heat at constant pressure 

are given by: 


t 

0 

c 

P 

200 F 

0.2255 Btu/ (lbm F) 

500 

0.1239 

800 

0.2093 

1100 

0.2091 

1400 

0.2120 


Po lynomia 1 Fit: 

Temperature Range: 660 to 1860 R 

Equation : 


c° = (0.2255 - 0.016292 0 + 0.00539 0 2 - 0.000758 6 3 + 
P 

0.000054 B 4 ) X 32.174 Btu/ (slug R) 


(B .44) 
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where 


T - 659.67 R 
300 R 


(B.44a) 


The variation of with pressure is given by 


= c - T 
P P 


/ 


3 el 


3T 


d P 


2 Vp 


Integrating along the isotherm T from p - 0 psig to p’ 
using Eqs. B.36 and B.38 results in 


= P: 


2T$ 2 

C P = c p “ TIT ln (1 + 


(B.45) 


Maximum Error: At zero pressure, the maximum error was 0.075%. 


6. Enthalpy 

The variation of enthalpy with both pressure and temperature is 
given by 

dh = c dT + -7 [l - TB] dp 
P P 


The enthalpy was arbitrarily chosen to be zero near the melting point, 
or T = 469.67°R. The above expression was integrated along an isobar 
p ss 0 psig from T* = 469.67°R to T 1 = T, and then along an isotherm 
T from p' = 0 psig to p 1 = p, to give 

6 2 0 3 

h « [300 (0.2255 6 - 0.016292 y + 0.00539 — - 

q4 q5 w mO 

0.000758 y + 0.000054 — ) + In (1 + xp) x (B.46) 

1 


778.26 


] 32.174 Btu/ slug 
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3 2 

where P and P are in lbm/ft , < is in ft /Ibf T in R, ® is given in 

° -1 2 
Eq. B.44a, 8 is in R and p is gage pressure in lbf/ft . 

7. Thermal Conductivity 
Data Points: 


t 

k 

200 F 

13.36 Btu/ (hr ft F) 

500 

14.57 

800 

15.18 

1100 

15.03 

1400 

14.13 


Polynomial Fit: 

Temperature Range: 660 to 1860 R 

Equation: 

k = (13.36 + 1.414167 0 - 0.142083 0 2 - 0.069167 0 3 + 

(B.47) 

0.007083 0 ; Btu/(hr ft R) 


where 


T - 659.67 R 
300 R 


Maximum Error: 0.2% 

8. Temperature Variation of Thermal Conductivity : 

Eq. B.47 was differentiated with respect to temperature to yield 


1 dk _ 1 (1. 414167 - 0.2841668 - 0.2075018 2 + 0.0283326 3 1 ( B .4 8 ) 

k dT 300 (13.36 + 1.4141670 - 0.1420838 2 - 0.0691679 3 + O.OO7O830 4 ) R 
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9. Dynamic Viscosity 
Data Points: 


t 

V 

200 F 

1.1316 Ibm/ (ft hr) 

500 

0.746 

800 

0.534 

1100 

0.411 

1400 

0.340 


Polynomial Fit: 

Temperature Range: 660 to 1860 R 

Equation: 


y = (1.316 - 0.896667 6 + 0.419833 6 2 - 0.102833 6 3 + 
0.009667 6 4 )/32.174 slug/(ft hr) 


where 


T - 659. 67R 
300R 


(B.49) 


Maximum Error : 1 . 2 % 
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IV. FC-75 INERT FLUOROCHEMICAL LIQUID 

The thermodynamic and transport properties of FC-75 fluid are ex- 
tracted from "3M Brand Inert Fluorochemical Liquids, 3M Company, Chemical 
Division, 1965." 

Due to the lack of experimental data, all properties were evaluated 
at atmospheric pressure and were assumed to be pressure independent. 

At one atmosphere some typical properties are: 

Nominal Boiling Point: 216 F 

Pour Point: -135 F 

Surface Tension, at 77F: 15 dynes/cm 

1, Equation of State 
Data Points: 


t 

P 

— 

- 50 F 

120.7 

Ibm/f t 3 

70 

110.5 


190 

100.3 



Polynomial Fit: 

Temperature Range: -80 to 216 F 

Equation : 

P = (155.522 - 0.085 T R~ 1 )x32.174 slug/ft 3 (B.50) 

Maximum Error: There was no difference between the computed and 

input data, within the accuracy of computation. 
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2. Isobaric Thermal Expansion Coefficient 

Using the definition of the isobaric thermal expansion coefficient, 
Eq. B.4, and Eq. B.50, one obtains 

0.085 1 

B = r t7 (B. 51) 

155.522 - 0.085 T R 


3. Isothermal Compressibility 

Since the equation of state (Eq. B.50) was assumed to be pressure 
independent, the isothermal compressibility defined by Eq. B.7 was 
assigned the value of zero 


k e 0 


4. Specific Heat at Constant Pressure 

Data Points: The variation of the zero pressure specific heat at con- 

stant pressure is given by: 


t 

c° 

P 

80° F 

0.2464 Btu/(lbm F) 

140 

0.2610 

200 

0.2756 


Polynomial Fit: 

Temperature Range: 70 to 210 F 

Equation: 

c° = (0.115082 + 2 .4333 x 10" 4 T R _1 ) x 32.174 
P 


Btu/ (slug R) 


(B. 52) 




142 


The variation of c with pressure is given by Eq. B.27. 
With p independent of pressure, Eq. B.27 was integrated along 
an isotherm T from p 1 = 0 psig to p 1 = p resulting in: 


c = 
P 



2 

23 T 


(B. 53) 


A check on the magnitude of the terms in Eq. B.53, using 

2 

typical running conditions , showed that the term 23 TP/p was 

only 0.000647o of c . Therefore, c was taken to be a function 

P P 

of temperature alone, namely 


c = c° (B. 54) 

P P 


Maximum Error: 0.025% 

5 . Enthalpy 

The variation of enthalpy with both pressure and temperature is 
given by 


dh = c° dT + - (1 - T3) dp (B.55) 

P P 

This expression was integrated along the isobar p = 0 psig from 
T* = T = 324.67 R to T' = T, and then along an isotherm T from p' = 

0 psig to p' = p to give 


h = [0.115082 (T - T ) + ~ x 2.4333 x 10~ 4 
L o 2 


(I 2 - T 2 )] x 32.174 + (1 - T 6) Btu/slug 


(B. 56) 


3 2 

where P is in slug/ft , p is gage pressure in Ibf /f t , T is in R and 3 
is in R \ 
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A check of the magnitude of the terms in Eq. B.56, using typical 
running conditions, showed that the last term which accounts for the 
pressure variation is only 0.00087% of h. Hence the enthalpy was taken 
to be a function of temperature alone, or 


h = [0.115082 (T - T ) + j x 2.4333 x 10“ 4 


(T 3 - T 3 )] x 32.174 Btu/slug 


(B. 57) 


Maximum Error: No data for enthalpy at atmospheric pressure were 
available for comparison. However, when the value of c^ 

(Eq. B.52) was compared with the result of differentiation of 
h with respect to temperature, there was no difference within 
the accuracy of the computation. 


6, Thermal Conductivity 
Data Points: 


t 

k 

- 50 F 

0.08745 Btu/ (hr ft F) 

50 

0.0809 

150 

0.0744 


Polynomial Fit: 

Temperature Range: -100 to 216 F 

Equation: 

k = 0.114181 - 6.53 x 10" 3 T r" 1 Btu/ (hr ft F) (B.58) 


Maximum Error: 0.0447, 
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7. Dynamic Viscosity 

Data Points: The variation of kinematic viscosity with temperature 

is given by: 



Polynomial Fit: 

Temperature Range : - 80 to 190 F 

Equation: 

v = e (1 „639 - 1.312933 6 + 0.25265 0 2 

3 

- 0.02471667 0 ) centistokes 


(B.59) 


where 


T - 409.67 R 


The dynamic viscosity is given by 


U » vp 


(B. 60) 


Maximum Error: 0.82% 
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V. FC-43 INERT FLUROCHEMICAL LIQUID 

The thermodynamic and transport properties of the FC-43 are ex- 
tracted from ”3M Brand Inert Fluorochemical Liquids, 3M Company, Chemical 
Division, 1965.” 

Due to the lack of experimental data, all properties were evalua- 
ted at atmospheric pressure and were assumed to be pressure independent. 

At one atmosphere some typical properties are: 

Nominal Boiling Point: 345 F 

Pour Point: - 58 F 

Surface Tension at 77F: 16 dynes/cm 

1. Equation of State 
Data Points : 


t 

P 


- 20 F 

123.6 

lbm/f t^ 

130 

112.15 


280 

100.75 



Polynomial Fit: 

Temperature Range : - 50 to 340 F 

Equation: 

p = (157.0883 - 0.076167 T R ) x 32.174 slug/ft (B.61) 

Maximum Error: There was no difference between the computed 

and input data, within the accuracy of computation. 
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2. Isobaric Thermal Expansion Coefficient 

Using the definition of the isobaric thermal expansion coefficient, 
Eq. B.4, and Eq. B.61, one obtains 


0.076167 

157.0883 - 0.076167 T R 


1 

R 


(B. 62) 


3, Isothermal Compressibility 

Since the equation of state (Eq. B.61) was assumed to be pressure 
independent, the isothermal compressibility defined by Eq. B.7 was 
assigned the value of zero. 


K = 0 


4. Specific Heat at Constant Pressure 

Data Points: The variation of zero pressure specific heat at 

constant pressure is given by: 



Polynomial Fit: 

Temperature Range: 40 to 77 F 

Equation : 

c° = (-0.020092 + 5.4054 x 10" 4 T R _1 ) x 32.174 


Btu/ (slug R) 


(B. 63) 
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A similar procedure to that used in the case of FC-75 (See 
Section B IV 4) has shown that to a close approximation c^ 
is pressure independent, or 


c 

P 


(B.64) 


Maximum Error: There was no difference between the computed and 

input data within the accuracy of computation. 

5. Enthalpy 

The variation of enthalpy with both pressure and temperature is 
given by Eq. B.55. Integration of this equation in a procedure similar 
to that followed for the coolant fluid FC-75 yields 


h = [-0.020092 (T - T ) + “ x 5.4054 x 10" 4 (T 2 - T 2 )] x 
L o 2 o J 

32.174 + — (1 - TS ) Btu/slug (B.65) 

//O.ZOp 

where T q = 401.67 R, p is in slug/ft^, T is in R, 8 is in R 1 , and 
p is gage pressure in lbf/ft 2 . 

A check of the magnitude of the terms in Eq. B.64, using 
typical running conditions, showed that the last term which accounts 
for the pressure variation is only 0.00084% of h. Therefore the enthalpy 
was taken to be a function of temperature alone or 

h = [-0.020092 (T - T ) + \ x 5.4054 x 10~ 4 (T 2 - T 2 ) ] x 
u o 2 o J 


32.174 Btu/slug 


(B.66) 
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Maximum Error: No data for enthalpy at atmospheric pressure were 

available for comparison. However, when the value of c 

P 

(Eq. B.63) was compared with the result of differentiation of 
h with respect to temperature, there was no difference within 
the accuracy of the computations. 

6. Thermal Conductivity 
Data Points: 


t 

k 

- 50 F 

0.0512 Btu/ (hr ft F) 

50 

0.0487 

150 

0.0462 


Polynomial Fit: 

Temperature Range : - 58 to 250 F 

Equation: 


k = 0.061442 - 2.5 x 10“ 5 T R _1 Btu/ (hr ft F) (B.67) 

Maximum Error: 0.11% 


7. Dynamic Viscosity 

Data Points: The variation of kinematic viscosity with tempera- 

ture is given by: 


t 

V 


r 

- 20 F 

15.80 

centistokes 

70 

2.84 


160 

0.855 


250 

0.35 
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Polynomial Fit: 

Temperature Range : - 80 to 320 F 

Equation: 


v = e (2.76 - 2.043483 0 + 0.362 6' 


- 0.034717 


e 3 ) 


centistokes 


(B.68) 


where 


e = 


T - 439.67 R 
90 R 


The dynamic viscosity is given by 


y = v p 


(B.69) 


Maximum Error : 1.0% 
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APPENDIX C 


Optical Properties 


Three optical properties are required in the radiative analysis 
discussed in Chapter 6, namely the total hemispherical emittance 


e(T) = 


V T) 


oc 

/ 


e x (T) E A (T)dX 


(C.l) 


and the two auxiliary functions (see Eqs. 6.10 and 6.11) 


XX(T, T„) = 


1, 2' E b (T 2 ) 


00 

/ •»«! 


) s x (T 2 ) E b ^ x (T 2 )dX 


(C.2) 


XXX(T,,T-,T,) = 


1 ’ 2 ’ 3 ' E b ( T 3 > 


00 

/ £a(t1 


) s x ( t 2 ) Ex (T 3 ) E b>x (T 3 ) dX (C.3) 


In view of the temperature independence of the spectral emittance 
for dielectrics, the two functions XX and XXX are functions of a single 
temperature, the temperature of the surface element represented by the last 
subscript on the lef-hand sides of Eqs. 6.10 and 6.11: 


XX (T) 



(T)dA 


(C .4) 


XXX(T) = 


00 




(T)dX 


(C.5) 
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I. SURFACE COATING Z-93 


The functions C.l, C.4 and C.5 of the previous section are evaluated 
for the zinc oxide/potassium silicate coating Z-93 on the basis of spectral 
reflectance data measured by IITRI and published in the NASA Contractor 
Report No. 1420, titled Emissivity Coatings for Low- Temperature Space 
Radiators, by G. R. Cunning ton, J. R. Grammer, and F. J. Smith, Lockheed 
Aircraft Corp., Sunnyvale, Calif., Sept. 1969, pp. 66 through 81. 

The evaluated functions defined through Eqs. C.l, 4 and 5 are 
collocated by power polynomials of this form 


f (T) = 



(C.6) 


For the total hemispherical emittance, a fourth degree power polynomial was 
found to be satisfactory with 

a Q = 0.8990103 

a ± = -0.1400633 x 10' 3 
a 2 = 0.387900 x 10" 6 

a 3 = -0.3937509 x 10” 9 
a, = 0.1015627 x 10" 12 


For the auxiliary functions XX and XXX the coefficients are 


XX 


XXX 


a rt = 

0.7804112 



0.6538383 



0 

a, = 

-0.5527205 

X 

10‘ 4 

0.1144374 

X 

10 

1 

a n = 

0.2530228 

X 

io ” 6 

-0.2432286 

X 

10 ' 

2 

a = 

-0.3229181 

X 

10" 9 

-0.1437500 

X 

10' 

a 4 = 

0.8854202 

X 

10" 13 

0.4947915 

X 

10 


-3 

-7 

-9 

-13 
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APPENDIX D 

I. The Fin- To- Tube Shape Factor 


A closed-form integration for the view factor of the fin with respect 
to the tube was carried out by Mr. Yao. This view factor occurs in Eqs. 6.15 
and 6.16. Only the final results are given here. 

The reader should recognize that some of the symbols defined below 
(Eqs . D.l through 6) apply only here. 

Let (x f , y f , z f ) designate the position of the center of an area 

element A on the fin and z the same on the tube. Let r , s and s represent, 
f m e l x 

respectively, the outer tube radius, the fin tip and the fin root thickness, 
and let the fin height be given as H. 

Then, with 

P ~ + y^ (D.l) 

f J f 


S = arc tg 



(D.2) 


a 


arc 



(D.3) 


0 - arc 



(D.4) 


2 2 

a = p + r + (z 


f “ Z m ) 


(D.5) 


$ 


arc 


r 

e 


P 


(D.6) 
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one obtains first 

a-2r P cos a-2r 2 cos (a+3) a-2r 2 cos (a+6) 

e e e 

^1 a-2r p cos(cf> -8) + a-2r p cos(<f> w -8) a-2r P cos (4>— B ) 


Z 2 = 


. 21 , ,2 ,,22 3 

8r P (z -z ) + 4P r a-a 

e t m e 

2(a 2 -4r 2 p 2 ) 3/2 

e 


(D.8) 


2r P -a cos(^*”8) 2r P -a cos ($-8) 

0 ? 6 

= arc sin — ~ arc sm — ~ rr — ttt 

3 a-2r P cos(<{> -8) a-2r P cos(<J>-8) 


Z 4 = 


2 2 2 2 
2a(z.-z )' - a Z + 4r P 
v f m e 

2 . 2 n 2 
a - 4r P 
e 


Z„ = 


r P sin(<J>*-8) r P sin (<l>-3) 
e e 


i-2r P cos (<J>*-8) a-2r P cos ($-8) 


The final result is 


-Z 


ss = aa A z 


1 sin(°H-ft) 


4Trp 


2irp 


2 + Z 2 Z 3 Z 4 Z 5 


(D.9) 


(D. 10) 


(D.ll) 


(D.12) 


This expression contains all the geometric relations that are required for 
fin-channel radiative interaction. It needs to be evaluated only once for 
every fin element. 

The values of SS given by Eq. D-12 vary greatly depending on the fin 
and tube elements under consideration. Typically the shape factor between 
tube element and adjacent fin element is three to five orders of magnitudes 
larger than the shape factor between tube element and the next closest fin 
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element. If this large variation in shape factor is allowed to remain, 
unrealistic oscillations in the fin radiosity will occur and sizeable truncation 
errors will result when the radiant fluxes are integrated across the fin 
surface . 

To eliminate these truncation errors the local shape factor for the 
root fin element is replaced by a mean valve between two adjacent elements 
which have the same Z location. The mean value for the shape factor between 
a tube element and its adjacent fin element is calculated from 


i-i . -ir- sin 20 -f 

ss -[i!!. y ] [. o+ — 5-2] 


(D. 13 ) t 


rather than Eq. D.12. 


Expressions for y and 0 q are: 


y = tan 


-1 


Zm - Z\ 
X, 


[cos <f)^ + tan (p sin - y 


r + 
e £ 


- tan 


o 
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II. Tube to Tube Shape Factor 


The simplified analysis has indicated that optimum dimensions of 
the radiator system will result in close tube spacing. Therefore for the 
case of an optiminually designed system, the radiant interaction between 
adjacent tubes must be taken into account. This section summarizes the 
results of the shape factor between two fin elements. The procedure used 
is Hottel’s crossed string method which is valid for infinitely long elements 
that are generated by a straight line moving parallel to itself. The finite 
length of the tube elements to accounted for by multiplying Hottel's result 
by a weighting factor. 

From Hottel's crossed string method (see Fig. 5 ) the shape factor 
between two infinitely long tubes is : 


? l-2 “["/2 1 - 0]) Y + [(4) 


2 2L ll/2 

+ COS 0) - 1 


- (IJ - cos 0 ) 


(D.14) 


where 




t = one half the fin thickness at its root. 

- outside radius of the tube 

L = one-half the distance between tube centers 


Y = 2 - < 6 ! + 0 1> 


0_ = tan 


-1 


2L - R 2 cos 0 


0, = COS 


-1 / R 2 sin 0 1 


The weighting factor for two tube elements of width Az located at 
Z^, Z 2 (see Fig. 5 ) is 
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6. 


■/ 


WF = j cos $d$ 

1 6 . 


\ <B 2 -6 1 ) +| (sin 2S 2 - sin 28 ^ 


(D. 15) 


where 


6 


1 


tan 


-1 


Z 2 - a Z/ 2 - Z x 
2(L-R 2 ) 


= tan 


-1 


Z 2 + AZ/2 - Z 1 
2(L-R^j 


The shape factor between tube element is now 


SS 


1-2 


WF 

tt/2 


1-2 


(D. 16) 


where the valve for F^ is given by Eq. D.14. 
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